NS Hemisphere Analysis
KE Lotterhos
10/27/2020
Based on basic code for calculating and mapping climatic novelty using the sigma dissimilarity metric. * written by Colin Mahony * Supp Material of “A closer look at novel climates: new methods and insights at continental to landscape scales”
What this code does
Compare 1800 analog (A) to today surface (B), with today’s surface variability (C)
Compare today surface analog (A) to 2100-RCP8.5 (B), with today’s surface variability (C)
Compare today surface analog (A) to 2100-RCP4.5 (B), with today’s surface variability (C)
This code also compares the converse (e.g. the current climate as B to the future climate as A) to assess disappearing climates
This code uses the annnual climate normals and monthly ICV for the novelty calculation
##################################
#### System setup ####
##################################
## Specify location of data ####
setwd("/Users/lotterhos/Documents/GitHub/OceanClimateNovelty")
#setwd("~/Desktop/PostDoc/OceanClimateNovelty/")
source("src/0-Novelty_Oceans_Functions.R")## Warning: package 'raster' was built under R version 3.6.2## Loading required package: sp## Warning: package 'sp' was built under R version 3.6.2## Loading required package: ade4## Loading required package: adehabitatMA## Registered S3 methods overwritten by 'adehabitatMA':
## method from
## print.SpatialPixelsDataFrame sp
## print.SpatialPixels sp##
## Attaching package: 'adehabitatMA'## The following object is masked from 'package:raster':
##
## buffer## Loading required package: CircStats## Loading required package: MASS##
## Attaching package: 'MASS'## The following objects are masked from 'package:raster':
##
## area, select## Loading required package: boot##
## Attaching package: 'data.table'## The following object is masked from 'package:raster':
##
## shift## ── Attaching packages ────────────────────── tidyverse 1.2.1 ──## ✓ ggplot2 3.3.2 ✓ purrr 0.3.3
## ✓ tibble 2.1.3 ✓ dplyr 0.8.3
## ✓ tidyr 1.0.0 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.4.0## Warning: package 'ggplot2' was built under R version 3.6.2## ── Conflicts ───────────────────────── tidyverse_conflicts() ──
## x dplyr::between() masks data.table::between()
## x tidyr::extract() masks raster::extract()
## x dplyr::filter() masks stats::filter()
## x dplyr::first() masks data.table::first()
## x dplyr::id() masks adehabitatLT::id()
## x dplyr::lag() masks stats::lag()
## x dplyr::last() masks data.table::last()
## x dplyr::select() masks MASS::select(), raster::select()
## x purrr::transpose() masks data.table::transpose()## Loading required package: spam## Loading required package: dotCall64## Loading required package: grid## Spam version 2.4-0 (2019-11-01) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.##
## Attaching package: 'spam'## The following objects are masked from 'package:base':
##
## backsolve, forwardsolve## Loading required package: maps##
## Attaching package: 'maps'## The following object is masked from 'package:purrr':
##
## map## See https://github.com/NCAR/Fields for
## an extensive vignette, other supplements and source code## rgdal: version: 1.4-7, (SVN revision 845)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.4.2, released 2019/06/28
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/gdal
## GDAL binary built with GEOS: FALSE
## Loaded PROJ.4 runtime: Rel. 5.2.0, September 15th, 2018, [PJ_VERSION: 520]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/proj
## Linking to sp version: 1.3-1## Warning: replacing previous import 'sf::st_make_valid' by
## 'lwgeom::st_make_valid' when loading 'tmap'## Registered S3 method overwritten by 'xts':
## method from
## as.zoo.xts zoo## Checking rgeos availability: TRUE## Warning: package 'sf' was built under R version 3.6.2## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0##
## Attaching package: 'fasterize'## The following object is masked from 'package:graphics':
##
## plot##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combine##
## Attaching package: 'reshape2'## The following object is masked from 'package:tidyr':
##
## smiths## The following objects are masked from 'package:data.table':
##
## dcast, melt##################################
#### Read in the input data ####
##################################
path <- "/Users/lotterhos/Google Drive/_LotterhosLab2021/LotLab_projects/0-grants_or_active/2018-Unfunded-OceanClimateNovelty/Data/Lotterhos/"
path2<- "/Users/lotterhos/Google Drive/_LotterhosLab2021/LotLab_projects/0-grants_or_active/2018-Unfunded-OceanClimateNovelty/Data/"
dat1 <- fread(paste0(path,"Katie_Temp_Arag_1800_2000.txt"), sep = ",")
dat4.5 <- fread(paste0(path,"Katie_Temp_Arag_2070_2100_RCP45.txt"), sep = ",")
dat8.5 <- fread(paste0(path,"Katie_Temp_Arag_2070_2100_RCP85.txt"), sep = ",")
ICV <- fread(paste0(path2, "ICV_Temp_Arag_pH_1960_2020.txt"))
results_path <- "/Users/lotterhos/Google Drive/_LotterhosLab2021/LotLab_projects/0-grants_or_active/2018-Unfunded-OceanClimateNovelty/Results/"
results_dir <- paste0(results_path, "HemisphereRestricted/")
head(dat1)## No Lon Lat Year Month SST Arag pH
## 1: 1 20.5 -69.5 1800 1 -0.05589 1.8167 8.0967
## 2: 1 20.5 -69.5 1800 2 -0.49722 2.0572 8.1706
## 3: 1 20.5 -69.5 1800 3 -0.88001 1.9321 8.1444
## 4: 1 20.5 -69.5 1800 4 -1.24800 1.8865 8.1396
## 5: 1 20.5 -69.5 1800 5 -2.23610 1.8170 8.1375
## 6: 1 20.5 -69.5 1800 6 -0.21900 1.9615 8.1438head(dat8.5)## No Lon Lat Year Month SST Arag pH
## 1: 1 20.5 -69.5 2070 1 0.503480 0.95276 7.7819
## 2: 1 20.5 -69.5 2070 2 0.062159 1.00710 7.8161
## 3: 1 20.5 -69.5 2070 3 -0.320630 0.96601 7.8035
## 4: 1 20.5 -69.5 2070 4 -0.688590 0.94601 7.8008
## 5: 1 20.5 -69.5 2070 5 -1.676700 0.90633 7.7981
## 6: 1 20.5 -69.5 2070 6 0.340380 0.98976 7.8048head(dat4.5)## No Lon Lat Year Month SST Arag pH
## 1: 1 20.5 -69.5 2070 1 0.58164 1.1518 7.8685
## 2: 1 20.5 -69.5 2070 2 0.14032 1.2362 7.9107
## 3: 1 20.5 -69.5 2070 3 -0.24248 1.1804 7.8955
## 4: 1 20.5 -69.5 2070 4 -0.61043 1.1555 7.8924
## 5: 1 20.5 -69.5 2070 5 -1.59860 1.1084 7.8899
## 6: 1 20.5 -69.5 2070 6 0.41854 1.2072 7.8965head(ICV)## No Lon Lat Year Month SST Arag pH
## 1: 1 20.5 -69.5 1960 1 -0.20420 1.5702 8.0291
## 2: 1 20.5 -69.5 1960 2 -0.64552 1.7479 8.0921
## 3: 1 20.5 -69.5 1960 3 -1.02830 1.6505 8.0698
## 4: 1 20.5 -69.5 1960 4 -1.39630 1.6131 8.0657
## 5: 1 20.5 -69.5 1960 5 -2.38440 1.5515 8.0635
## 6: 1 20.5 -69.5 1960 6 -0.36730 1.6803 8.0699unique(dat1$Year)## [1] 1800 1810 1820 1830 1970 1980 1990 2000unique(dat8.5$Year)## [1] 2070 2080 2090 2100unique(dat4.5$Year)## [1] 2070 2080 2090 2100unique(ICV$Year)## [1] 1960 1970 1980 1990 2000 2010 2020hist(c(dat1$Arag, dat8.5$Arag))
hist(log10(c(dat1$Arag, dat8.5$Arag)))
boxplot ICV
pdf(paste0(results_dir, "ICV_Lat.pdf"), width=9, height=7)
par(mfrow=c(3,1), mar=c(4,5,0.5,0.5), oma=c(2,0,0,0), las=0)
lat <- sort(unique(ICV$Lat))
length(lat)## [1] 169lat[1]## [1] -78.5lat[42]## [1] -37.5lat[80]## [1] 0.5lat[117]## [1] 37.5lat[158]## [1] 78.5ind <- c(1,42,80,117,158)
lat_names=rep(NA, 169)
lat_names[ind] <- lat[ind]
lat_names[80] <- 0 #close enough
boxplot(ICV$SST~ICV$Lat, col="green",
#border = rgb(0,0,0,0.3), outpch=20, outcol= rgb(0,0,0,0.3), outbg=rgb(0,0,0,0.3),
ylab="", xlab="", names=lat_names, las=1, cex=1.2)
mtext("SST", side=2, cex=1.5, line=3)
boxplot(ICV$Arag~ICV$Lat, col="green",
ylab="", xlab="", names=lat_names, las=1, cex=1.2)
mtext("Arag.", side=2, cex=1.5, line=3)
boxplot(ICV$pH ~ ICV$Lat, col="green",
ylab="", xlab="", names=lat_names, las=1, cex=1.2)
mtext("pH", side=2, cex=1.5, line=3)
mtext("Latitude", side=1, outer=TRUE, cex=2)
dev.off()## quartz_off_screen
## 2##################################
#### Log-transform Arag ####
##################################
dat1$Arag <- log10(dat1$Arag)
dat4.5$Arag <- log10(dat4.5$Arag)
dat8.5$Arag <- log10(dat8.5$Arag)
ICV$Arag <- log10(ICV$Arag)
# Aragonite saturation is a ratio variable; it is limited at zero and
# proportional changes are meaningful. The analysis uses raw values of
# aragonite saturation state. Under raw scaling, the difference between 1 and 2
# is equal to the difference between zero and one, which doesn’t reflect the
# meaning of the variable. Instead, the difference between 1 and 2 should
# be given the same significance as the difference between 0.5 and 1 (doubling vs halving).
# Log-transforming aragonite saturation state will produce this more meaningful
# proportional scaling.
# More generally on the same topic: The need for proportional scaling is the
# reason why precipitation and other ratio variables (e.g. degree-days) are
# typically log-transformed in climate space analysis. In theory, all variables
# should be log-transformed to produce an environmental space where distances
# are comparable for different locations; however, in practice temperature
# doesn’t need to be log-transformed because it doesn’t vary across orders
# of magnitude, and in our case pH is already a log-scaled variable.
#boxplot(dat8.5$pH ~ dat8.5$Year)
#boxplot(dat4.5$pH ~ dat4.5$Year)
# note that these two datasets are exactly the same for years < 2000#--------------------------------
#### 40-year climate normals ####
#--------------------------------
# for example, 1930 represents from 1/1/1925 to 12/31/1934.
dat_2000 <- dat1 %>% filter(Year>1960 & Year<2010)
dim(dat_2000)## [1] 1903776 8# Below is a standardization code that can be deleted becuase the PCA
# does the standardization
# x_SST <- mean(dat_2000$SST)
# s_SST <- sd(dat_2000$SST)
# x_Arag <- mean(dat_2000$Arag)
# s_Arag <- sd(dat_2000$Arag)
# x_pH <- mean(dat_2000$pH)
# s_pH <- sd(dat_2000$pH)
#
# head(dat_2000)
# dat_2000$SST <- (dat_2000$SST - x_SST)/s_SST
# dat_2000$Arag <- (dat_2000$Arag - x_Arag)/s_Arag
# dat_2000$pH <- (dat_2000$pH - x_pH)/s_pH
# dat_2000 %>% summarise_each(sd)
#
dat_1800 <- dat1 %>% filter(Year<1850)
dim(dat_1800)## [1] 1903776 8# dat_1800$SST <- (dat_1800$SST - x_SST)/s_SST
# dat_1800$Arag <- (dat_1800$Arag - x_Arag)/s_Arag
# dat_1800$pH <- (dat_1800$pH - x_pH)/s_pH
#
#
dat_2100_8.5 <- dat8.5
# dim(dat_2100_8.5)
# dat_2100_8.5$SST <- (dat_2100_8.5$SST - x_SST)/s_SST
# dat_2100_8.5$Arag <- (dat_2100_8.5$Arag - x_Arag)/s_Arag
# dat_2100_8.5$pH <- (dat_2100_8.5$pH - x_pH)/s_pH
#
dat_2100_4.5 <- dat4.5
# dim(dat_2100_4.5)
# dat_2100_4.5$SST <- (dat_2100_4.5$SST - x_SST)/s_SST
# dat_2100_4.5$Arag <- (dat_2100_4.5$Arag - x_Arag)/s_Arag
# dat_2100_4.5$pH <- (dat_2100_4.5$pH - x_pH)/s_pH
sum(!complete.cases(dat_1800))## [1] 0sum(!complete.cases(dat_2000))## [1] 0sum(!complete.cases(dat_2100_4.5))## [1] 0sum(!complete.cases(dat_2100_8.5))## [1] 0sum(!complete.cases(ICV))## [1] 0# Calculate climate normals
# Returns the mean climate for each station at each time
norm_1800 <- calculate_normals_annual(dat_1800)
norm_2000 <- calculate_normals_annual(dat_2000)
norm_2100_8.5 <- calculate_normals_annual(dat_2100_8.5)
norm_2100_4.5 <- calculate_normals_annual(dat_2100_4.5)
head(norm_1800)## No SST Arag pH
## 1 1 -1.2720933 0.2661633 8.129290
## 2 2 -1.1802993 0.2581453 8.136665
## 3 3 -1.0816947 0.2502643 8.141752
## 4 4 -0.9757862 0.2425956 8.144575
## 5 5 -0.8618450 0.2352653 8.145356
## 6 6 -0.7388042 0.2284502 8.144419dim(norm_1800)## [1] 39662 4dim(norm_2000)## [1] 39662 4dim(norm_2100_8.5)## [1] 39662 4dim(norm_2100_4.5)## [1] 39662 4#--------------------------------
### data frame to link station number to Lat Long ####
#--------------------------------
head(dat8.5)## No Lon Lat Year Month SST Arag pH
## 1: 1 20.5 -69.5 2070 1 0.503480 -0.021016484 7.7819
## 2: 1 20.5 -69.5 2070 2 0.062159 0.003072596 7.8161
## 3: 1 20.5 -69.5 2070 3 -0.320630 -0.015018378 7.8035
## 4: 1 20.5 -69.5 2070 4 -0.688590 -0.024104273 7.8008
## 5: 1 20.5 -69.5 2070 5 -1.676700 -0.042713644 7.7981
## 6: 1 20.5 -69.5 2070 6 0.340380 -0.004470102 7.8048stations <- unique(dat8.5$No)
stationInfo = data.frame(stations=stations, lat=NA, long=NA)
unik <- which(!duplicated(dat8.5$No))
head(unik)## [1] 1 49 97 145 193 241stationInfo$lat <- dat8.5$Lat[unik]
stationInfo$long <- dat8.5$Lon[unik]
names(stationInfo)[1] <- "No"
head(stationInfo)## No lat long
## 1 1 -69.5 20.5
## 2 2 -68.5 20.5
## 3 3 -67.5 20.5
## 4 4 -66.5 20.5
## 5 5 -65.5 20.5
## 6 6 -64.5 20.5which(!complete.cases(stationInfo))## integer(0)#--------------------------------
### ICV matrix
head(ICV) ## No Lon Lat Year Month SST Arag pH
## 1: 1 20.5 -69.5 1960 1 -0.20420 0.1959550 8.0291
## 2: 1 20.5 -69.5 1960 2 -0.64552 0.2425166 8.0921
## 3: 1 20.5 -69.5 1960 3 -1.02830 0.2176155 8.0698
## 4: 1 20.5 -69.5 1960 4 -1.39630 0.2076613 8.0657
## 5: 1 20.5 -69.5 1960 5 -2.38440 0.1907518 8.0635
## 6: 1 20.5 -69.5 1960 6 -0.36730 0.2253868 8.0699C <- ICV %>% select(No, SST, Arag, pH)
head(C)## No SST Arag pH
## 1: 1 -0.20420 0.1959550 8.0291
## 2: 1 -0.64552 0.2425166 8.0921
## 3: 1 -1.02830 0.2176155 8.0698
## 4: 1 -1.39630 0.2076613 8.0657
## 5: 1 -2.38440 0.1907518 8.0635
## 6: 1 -0.36730 0.2253868 8.0699#-------------------------------- #---------------
### Hemisphere Indexes ###
# N. hemisphere
N_hem <- which(stationInfo$lat>=0)
S_hem <- which(stationInfo$lat<0)
#---------------1800 analog to today (Novelty)
#--------------------------------
### 1800 analog to today (Novelty) ####
# What are today's novel climates compared to 1800?
# Novel climates are identified by comparing the
# future climate realization (B) for each gridpoint
# to the past (A) climate realizations (Williams)
# RESTRICT BASELINE TO HEMISPHERE OF OBSERVATION
#--------------------------------
# A is baseline, the base that a station from B is compared to
# If A is past and B is future, then degree of novelty
# If A is future and B is past, then degree of disappearance
# B will be subset to a particular station for analysis
# C is the data frame used to calculate the ICV, will also be subset to a particular station for analysis
A <- norm_1800
# 1800 climate normals
# Serves as the Baseline
head(A)## No SST Arag pH
## 1 1 -1.2720933 0.2661633 8.129290
## 2 2 -1.1802993 0.2581453 8.136665
## 3 3 -1.0816947 0.2502643 8.141752
## 4 4 -0.9757862 0.2425956 8.144575
## 5 5 -0.8618450 0.2352653 8.145356
## 6 6 -0.7388042 0.2284502 8.144419dim(A)## [1] 39662 4# subset the 2000 data to the same stations as the 1800 data
B <- norm_2000
# 2000 climate normals
head(B)## No SST Arag pH
## 1 1 -1.4479254 0.1799026 8.035921
## 2 2 -1.3696480 0.1741309 8.045892
## 3 3 -1.2845583 0.1683092 8.053425
## 4 4 -1.1916660 0.1624748 8.058488
## 5 5 -1.0895625 0.1567233 8.061231
## 6 6 -0.9764333 0.1512065 8.062000dim(B)## [1] 39662 4# same columns as A
# sanity check to make sure stations in right order
identical(A$No, B$No) # should be true## [1] TRUE# check list 18906,18952, 29736, 29789, 8193, 29319,
# calc_sigma_D(A, B, C, 18906, "_1800_today")
# calc_sigma_D(A, B, C, 18952, "_1800_today")
# calc_sigma_D(A, B, C, 29736, "_1800_today")
# calc_sigma_D(A, B, C, 29789, "_1800_today")
# calc_sigma_D(A, B, C, 8193, "_1800_today")
# calc_sigma_D(A, B, C, 29319, "_1800_today")
#NN.sigma_1800_today <- loop_sigma_D(A, B, C, "_1800_today") # uses whole ocean as a baseline
NN.sigma_1800_today_N <- loop_sigma_D(A[N_hem,], B[N_hem,], C, "_1800_today")## [1] 100 317 17413
## [1] 200 610 17413
## [1] 300 934 17413
## [1] 400 1247 17413
## [1] 500 1683 17413
## [1] 600 2053 17413
## [1] 700 2356 17413
## [1] 800 2660 17413
## [1] 900 2971 17413
## [1] 1000 3213 17413
## [1] 1100 3452 17413
## [1] 1200 3759 17413
## [1] 1300 3996 17413
## [1] 1400 4236 17413
## [1] 1500 4540 17413
## [1] 1600 4847 17413
## [1] 1700 5220 17413
## [1] 1800 5530 17413
## [1] 1900 5830 17413
## [1] 2000 6129 17413
## [1] 2100 6427 17413
## [1] 2200 6726 17413
## [1] 2300 7099 17413
## [1] 2400 7694 17413
## [1] 2500 8267 17413
## [1] 2600 8629 17413
## [1] 2700 8887 17413
## [1] 2800 9069 17413
## [1] 2900 9286 17413
## [1] 3000 9454 17413
## [1] 3100 9622 17413
## [1] 3200 9790 17413
## [1] 3300 9954 17413
## [1] 3400 10119 17413
## [1] 3500 10287 17413
## [1] 3600 10458 17413
## [1] 3700 10592 17413
## [1] 3800 10764 17413
## [1] 3900 10945 17413
## [1] 4000 11095 17413
## [1] 4100 11309 17413
## [1] 4200 11469 17413
## [1] 4300 11702 17413
## [1] 4400 11871 17413
## [1] 4500 12109 17413
## [1] 4600 12278 17413
## [1] 4700 12450 17413
## [1] 4800 12692 17413
## [1] 4900 12863 17413
## [1] 5000 13035 17413
## [1] 5100 13280 17413
## [1] 5200 13453 17413
## [1] 5300 13627 17413
## [1] 5400 13878 17413
## [1] 5500 14056 17413
## [1] 5600 14235 17413
## [1] 5700 14414 17413
## [1] 5800 14670 17413
## [1] 5900 14848 17413
## [1] 6000 15026 17413
## [1] 6100 15204 17413
## [1] 6200 15384 17413
## [1] 6300 15641 17413
## [1] 6400 15819 17413
## [1] 6500 15998 17413
## [1] 6600 16176 17413
## [1] 6700 16356 17413
## [1] 6800 16534 17413
## [1] 6900 16791 17413
## [1] 7000 16970 17413
## [1] 7100 17149 17413
## [1] 7200 17328 17413
## [1] 7300 17506 17413
## [1] 7400 17684 17413
## [1] 7500 17862 17413
## [1] 7600 18119 17413
## [1] 7700 18297 17413
## [1] 7800 18476 17413
## [1] 7900 18655 17413## [1] 8000 18912 17413## [1] 8100 19090 17413
## [1] 8200 19269 17413
## [1] 8300 19526 17413
## [1] 8400 19704 17413
## [1] 8500 19881 17413
## [1] 8600 20058 17413
## [1] 8700 20312 17413
## [1] 8800 20488 17413
## [1] 8900 20664 17413
## [1] 9000 20916 17413
## [1] 9100 21091 17413
## [1] 9200 21266 17413
## [1] 9300 21441 17413
## [1] 9400 21691 17413
## [1] 9500 21866 17413
## [1] 9600 22041 17413
## [1] 9700 22218 17413
## [1] 9800 22468 17413
## [1] 9900 22643 17413
## [1] 10000 22817 17413
## [1] 10100 22992 17413
## [1] 10200 23240 17413
## [1] 10300 23414 17413
## [1] 10400 23662 17413
## [1] 10500 23836 17413
## [1] 10600 24009 17413
## [1] 10700 24256 17413
## [1] 10800 24430 17413
## [1] 10900 24679 17413
## [1] 11000 24927 17413
## [1] 11100 25250 17413
## [1] 11200 25500 17413
## [1] 11300 25823 17413
## [1] 11400 26075 17413
## [1] 11500 26397 17413
## [1] 11600 26790 17413
## [1] 11700 27106 17413
## [1] 11800 27494 17413
## [1] 11900 27959 17413
## [1] 12000 28351 17413
## [1] 12100 28892 17413
## [1] 12200 29303 17413## [1] 12300 29551 17413
## [1] 12400 29721 17413## [1] 12500 29857 17413
## [1] 12600 30024 17413
## [1] 12700 30174 17413
## [1] 12800 30333 17413
## [1] 12900 30505 17413
## [1] 13000 30643 17413
## [1] 13100 30827 17413
## [1] 13200 30971 17413
## [1] 13300 31118 17413
## [1] 13400 31311 17413
## [1] 13500 31458 17413
## [1] 13600 31658 17413
## [1] 13700 31811 17413
## [1] 13800 32021 17413
## [1] 13900 32177 17413
## [1] 14000 32392 17413
## [1] 14100 32552 17413
## [1] 14200 32775 17413
## [1] 14300 32941 17413
## [1] 14400 33112 17413
## [1] 14500 33355 17413
## [1] 14600 33534 17413
## [1] 14700 33790 17413
## [1] 14800 33968 17413
## [1] 14900 34146 17413
## [1] 15000 34400 17413
## [1] 15100 34577 17413
## [1] 15200 34753 17413
## [1] 15300 35004 17413
## [1] 15400 35179 17413
## [1] 15500 35353 17413
## [1] 15600 35602 17413
## [1] 15700 35777 17413
## [1] 15800 35952 17413
## [1] 15900 36129 17413
## [1] 16000 36377 17413
## [1] 16100 36553 17413
## [1] 16200 36728 17413
## [1] 16300 36975 17413
## [1] 16400 37148 17413
## [1] 16500 37391 17413
## [1] 16600 37635 17413
## [1] 16700 37877 17413
## [1] 16800 38188 17413
## [1] 16900 38498 17413
## [1] 17000 38808 17413
## [1] 17100 39119 17413
## [1] 17200 39428 17413
## [1] 17300 39759 17413
## [1] 17400 40107 17413NN.sigma_1800_today_S <- loop_sigma_D(A[S_hem,], B[S_hem,], C, "_1800_today")## [1] 100 140 22249
## [1] 200 298 22249
## [1] 300 459 22249
## [1] 400 622 22249
## [1] 500 765 22249
## [1] 600 909 22249
## [1] 700 1052 22249
## [1] 800 1195 22249
## [1] 900 1345 22249
## [1] 1000 1468 22249
## [1] 1100 1619 22249
## [1] 1200 1763 22249
## [1] 1300 1923 22249
## [1] 1400 2093 22249
## [1] 1500 2228 22249
## [1] 1600 2401 22249
## [1] 1700 2539 22249
## [1] 1800 2718 22249
## [1] 1900 2858 22249
## [1] 2000 3040 22249
## [1] 2100 3182 22249
## [1] 2200 3372 22249
## [1] 2300 3518 22249
## [1] 2400 3710 22249
## [1] 2500 3856 22249
## [1] 2600 4047 22249
## [1] 2700 4190 22249
## [1] 2800 4368 22249
## [1] 2900 4506 22249
## [1] 3000 4678 22249
## [1] 3100 4810 22249
## [1] 3200 4967 22249
## [1] 3300 5092 22249
## [1] 3400 5240 22249
## [1] 3500 5367 22249
## [1] 3600 5531 22249
## [1] 3700 5666 22249
## [1] 3800 5834 22249
## [1] 3900 5969 22249
## [1] 4000 6141 22249
## [1] 4100 6278 22249
## [1] 4200 6451 22249
## [1] 4300 6586 22249
## [1] 4400 6753 22249
## [1] 4500 6882 22249
## [1] 4600 7031 22249
## [1] 4700 7149 22249
## [1] 4800 7275 22249
## [1] 4900 7388 22249
## [1] 5000 7514 22249
## [1] 5100 7640 22249
## [1] 5200 7753 22249
## [1] 5300 7882 22249
## [1] 5400 8010 22249
## [1] 5500 8139 22249## [1] 5600 8256 22249
## [1] 5700 8404 22249
## [1] 5800 8538 22249
## [1] 5900 8680 22249
## [1] 6000 8864 22249
## [1] 6100 9073 22249
## [1] 6200 9259 22249
## [1] 6300 9501 22249
## [1] 6400 9758 22249
## [1] 6500 10019 22249
## [1] 6600 10336 22249
## [1] 6700 10614 22249
## [1] 6800 10899 22249
## [1] 6900 11128 22249
## [1] 7000 11363 22249
## [1] 7100 11532 22249
## [1] 7200 11769 22249
## [1] 7300 11938 22249
## [1] 7400 12178 22249
## [1] 7500 12350 22249
## [1] 7600 12523 22249
## [1] 7700 12772 22249
## [1] 7800 12948 22249
## [1] 7900 13203 22249
## [1] 8000 13382 22249
## [1] 8100 13561 22249
## [1] 8200 13824 22249
## [1] 8300 14012 22249
## [1] 8400 14279 22249
## [1] 8500 14462 22249
## [1] 8600 14649 22249
## [1] 8700 14915 22249
## [1] 8800 15102 22249
## [1] 8900 15289 22249
## [1] 9000 15556 22249
## [1] 9100 15740 22249
## [1] 9200 15926 22249
## [1] 9300 16112 22249
## [1] 9400 16386 22249
## [1] 9500 16573 22249
## [1] 9600 16760 22249
## [1] 9700 16947 22249
## [1] 9800 17224 22249
## [1] 9900 17414 22249
## [1] 10000 17604 22249
## [1] 10100 17881 22249
## [1] 10200 18065 22249
## [1] 10300 18247 22249
## [1] 10400 18428 22249
## [1] 10500 18688 22249
## [1] 10600 18867 22249
## [1] 10700 19044 22249
## [1] 10800 19298 22249
## [1] 10900 19475 22249
## [1] 11000 19652 22249
## [1] 11100 19828 22249
## [1] 11200 20081 22249
## [1] 11300 20258 22249
## [1] 11400 20436 22249
## [1] 11500 20694 22249
## [1] 11600 20874 22249
## [1] 11700 21054 22249
## [1] 11800 21314 22249
## [1] 11900 21494 22249
## [1] 12000 21674 22249
## [1] 12100 21934 22249
## [1] 12200 22115 22249
## [1] 12300 22295 22249
## [1] 12400 22553 22249
## [1] 12500 22731 22249
## [1] 12600 22909 22249
## [1] 12700 23161 22249
## [1] 12800 23335 22249
## [1] 12900 23508 22249
## [1] 13000 23753 22249
## [1] 13100 23923 22249
## [1] 13200 24092 22249
## [1] 13300 24326 22249
## [1] 13400 24485 22249
## [1] 13500 24690 22249
## [1] 13600 24839 22249
## [1] 13700 24987 22249
## [1] 13800 25178 22249
## [1] 13900 25322 22249
## [1] 14000 25465 22249
## [1] 14100 25646 22249
## [1] 14200 25785 22249
## [1] 14300 25923 22249
## [1] 14400 26095 22249
## [1] 14500 26228 22249
## [1] 14600 26360 22249
## [1] 14700 26522 22249
## [1] 14800 26652 22249
## [1] 14900 26811 22249
## [1] 15000 26938 22249
## [1] 15100 27065 22249
## [1] 15200 27218 22249
## [1] 15300 27344 22249
## [1] 15400 27470 22249
## [1] 15500 27619 22249
## [1] 15600 27743 22249
## [1] 15700 27890 22249
## [1] 15800 28011 22249
## [1] 15900 28132 22249
## [1] 16000 28274 22249
## [1] 16100 28394 22249
## [1] 16200 28514 22249
## [1] 16300 28648 22249
## [1] 16400 28764 22249
## [1] 16500 28898 22249
## [1] 16600 29022 22249
## [1] 16700 29172 22249
## [1] 16800 29327 22249
## [1] 16900 29498 22249
## [1] 17000 29708 22249
## [1] 17100 30060 22249
## [1] 17200 30380 22249
## [1] 17300 30680 22249
## [1] 17400 30923 22249
## [1] 17500 31234 22249
## [1] 17600 31474 22249
## [1] 17700 31716 22249
## [1] 17800 31954 22249
## [1] 17900 32192 22249
## [1] 18000 32359 22249
## [1] 18100 32598 22249
## [1] 18200 32844 22249
## [1] 18300 33019 22249
## [1] 18400 33266 22249
## [1] 18500 33439 22249
## [1] 18600 33616 22249
## [1] 18700 33863 22249
## [1] 18800 34038 22249
## [1] 18900 34213 22249
## [1] 19000 34463 22249
## [1] 19100 34638 22249
## [1] 19200 34814 22249
## [1] 19300 35066 22249
## [1] 19400 35242 22249
## [1] 19500 35419 22249
## [1] 19600 35678 22249
## [1] 19700 35859 22249
## [1] 19800 36042 22249
## [1] 19900 36312 22249
## [1] 20000 36490 22249
## [1] 20100 36664 22249
## [1] 20200 36908 22249
## [1] 20300 37078 22249
## [1] 20400 37310 22249
## [1] 20500 37464 22249
## [1] 20600 37612 22249
## [1] 20700 37801 22249
## [1] 20800 37943 22249
## [1] 20900 38122 22249
## [1] 21000 38258 22249
## [1] 21100 38428 22249
## [1] 21200 38563 22249
## [1] 21300 38697 22249
## [1] 21400 38864 22249
## [1] 21500 38997 22249
## [1] 21600 39161 22249
## [1] 21700 39292 22249
## [1] 21800 39446 22249
## [1] 21900 39572 22249
## [1] 22000 39721 22249
## [1] 22100 39866 22249
## [1] 22200 40030 22249NN.sigma_1800_today <- rbind(NN.sigma_1800_today_N, NN.sigma_1800_today_S)
head(NN.sigma_1800_today)## No NN.sigma_1800_today NN.station_1800_today NN.Mdist_1800_today
## 1 36 0.033725924 30751 0.23355016
## 2 37 0.015960753 30751 0.16010043
## 3 38 0.010016616 31111 0.12668130
## 4 39 0.001600598 14563 0.05055509
## 5 40 0.018324556 1958 0.17162744
## 6 41 0.004808260 17177 0.08767895
## numPCs_1800_today
## 1 2
## 2 2
## 3 2
## 4 2
## 5 2
## 6 2hist(NN.sigma_1800_today$numPCs_1800_today)
final_dat <- full_join(stationInfo, NN.sigma_1800_today, by="No")
head(final_dat)## No lat long NN.sigma_1800_today NN.station_1800_today NN.Mdist_1800_today
## 1 1 -69.5 20.5 0.005340623 40066 0.0924152
## 2 2 -68.5 20.5 0.007501424 40009 0.1095737
## 3 3 -67.5 20.5 0.009629160 40067 0.1241975
## 4 4 -66.5 20.5 0.015490226 40009 0.1577081
## 5 5 -65.5 20.5 0.024201364 40066 0.1974682
## 6 6 -64.5 20.5 0.049860481 40066 0.2848815
## numPCs_1800_today
## 1 2
## 2 2
## 3 2
## 4 2
## 5 2
## 6 2final_dat[which(is.infinite(final_dat$NN.sigma_1800_today)),]## [1] No lat long
## [4] NN.sigma_1800_today NN.station_1800_today NN.Mdist_1800_today
## [7] numPCs_1800_today
## <0 rows> (or 0-length row.names)final_dat[which(!complete.cases(final_dat)),]## [1] No lat long
## [4] NN.sigma_1800_today NN.station_1800_today NN.Mdist_1800_today
## [7] numPCs_1800_today
## <0 rows> (or 0-length row.names)#NN.sigma_1800_today[which(is.infinite(NN.sigma_1800_today))] <- NA
identical(nrow(final_dat), nrow(stationInfo)) ## [1] TRUE# should be true, same number as stations as stationInfo
# Visualize
world <- map_data("world2")
Plot_nonInt(final_dat$lat, final_dat$long,
final_dat$NN.sigma_1800_today, world, "sigD novel 2000")
# The NN.sigma_1800_today represents the novelty of today's
# climate compared to 1800
# NN.dist_1800_today represents the Mahalanobis distance from the climate
# in today to it's nearest neighbor in all 1800 data
# No--> NN.station_1800_today represents the station that was queried 2000 (B)
# and it's nearest data in the global reference NN.station_1800_today (1800, A)
# SEE CODE BELOW FOR ADDITION OF NN SITESToday analog to 1800 Disappearance
#--------------------------------
### Today analog to 1800 Disappearance####
# What are 1800 disappearing climates?
# disappearing climates are identified by comparing
# each past gridpoint (A) to all future climate realizations (B)
#--------------------------------
calc_sigma_D(A, B, C, 1)## NN.sigma NN.station NN.Mdist num_PCs
## 1 0.005073206 22827 0.09006697 2calc_sigma_D(B, A, C, 1)## NN.sigma NN.station NN.Mdist num_PCs
## 1 0.0005525579 38202 0.02969763 2B <- norm_1800
A <- norm_2000
identical(A$No, B$No) # should be true## [1] TRUEidentical(sort(A$No), A$No) # should be true## [1] TRUE# check list 18906,18952, 29736, 29789, 8193, 29319,
# calc_sigma_D(A, B, C, 18906, "_today_1800")
# calc_sigma_D(A, B, C, 18952, "_today_1800")
# calc_sigma_D(A, B, C, 29736, "_today_1800")
# calc_sigma_D(A, B, C, 29789, "_today_1800")
# calc_sigma_D(A, B, C, 8193, "_today_1800")
# calc_sigma_D(A, B, C, 29319, "_today_1800")
NN.sigma_today_1800_N <- loop_sigma_D(A[N_hem,], B[N_hem,], C, "_today_1800")## [1] 100 317 17413
## [1] 200 610 17413
## [1] 300 934 17413
## [1] 400 1247 17413
## [1] 500 1683 17413
## [1] 600 2053 17413
## [1] 700 2356 17413
## [1] 800 2660 17413
## [1] 900 2971 17413
## [1] 1000 3213 17413
## [1] 1100 3452 17413
## [1] 1200 3759 17413
## [1] 1300 3996 17413
## [1] 1400 4236 17413
## [1] 1500 4540 17413
## [1] 1600 4847 17413
## [1] 1700 5220 17413
## [1] 1800 5530 17413
## [1] 1900 5830 17413
## [1] 2000 6129 17413
## [1] 2100 6427 17413
## [1] 2200 6726 17413
## [1] 2300 7099 17413
## [1] 2400 7694 17413
## [1] 2500 8267 17413
## [1] 2600 8629 17413
## [1] 2700 8887 17413
## [1] 2800 9069 17413
## [1] 2900 9286 17413
## [1] 3000 9454 17413
## [1] 3100 9622 17413
## [1] 3200 9790 17413
## [1] 3300 9954 17413
## [1] 3400 10119 17413
## [1] 3500 10287 17413
## [1] 3600 10458 17413
## [1] 3700 10592 17413
## [1] 3800 10764 17413
## [1] 3900 10945 17413
## [1] 4000 11095 17413
## [1] 4100 11309 17413
## [1] 4200 11469 17413
## [1] 4300 11702 17413
## [1] 4400 11871 17413
## [1] 4500 12109 17413
## [1] 4600 12278 17413
## [1] 4700 12450 17413
## [1] 4800 12692 17413
## [1] 4900 12863 17413
## [1] 5000 13035 17413
## [1] 5100 13280 17413
## [1] 5200 13453 17413
## [1] 5300 13627 17413
## [1] 5400 13878 17413
## [1] 5500 14056 17413
## [1] 5600 14235 17413
## [1] 5700 14414 17413
## [1] 5800 14670 17413
## [1] 5900 14848 17413
## [1] 6000 15026 17413
## [1] 6100 15204 17413
## [1] 6200 15384 17413
## [1] 6300 15641 17413
## [1] 6400 15819 17413
## [1] 6500 15998 17413
## [1] 6600 16176 17413
## [1] 6700 16356 17413
## [1] 6800 16534 17413
## [1] 6900 16791 17413
## [1] 7000 16970 17413
## [1] 7100 17149 17413
## [1] 7200 17328 17413
## [1] 7300 17506 17413
## [1] 7400 17684 17413
## [1] 7500 17862 17413
## [1] 7600 18119 17413
## [1] 7700 18297 17413
## [1] 7800 18476 17413
## [1] 7900 18655 17413## [1] 8000 18912 17413## [1] 8100 19090 17413
## [1] 8200 19269 17413
## [1] 8300 19526 17413
## [1] 8400 19704 17413
## [1] 8500 19881 17413
## [1] 8600 20058 17413
## [1] 8700 20312 17413
## [1] 8800 20488 17413
## [1] 8900 20664 17413
## [1] 9000 20916 17413
## [1] 9100 21091 17413
## [1] 9200 21266 17413
## [1] 9300 21441 17413
## [1] 9400 21691 17413
## [1] 9500 21866 17413
## [1] 9600 22041 17413
## [1] 9700 22218 17413
## [1] 9800 22468 17413
## [1] 9900 22643 17413
## [1] 10000 22817 17413
## [1] 10100 22992 17413
## [1] 10200 23240 17413
## [1] 10300 23414 17413
## [1] 10400 23662 17413
## [1] 10500 23836 17413
## [1] 10600 24009 17413
## [1] 10700 24256 17413
## [1] 10800 24430 17413
## [1] 10900 24679 17413
## [1] 11000 24927 17413
## [1] 11100 25250 17413
## [1] 11200 25500 17413
## [1] 11300 25823 17413
## [1] 11400 26075 17413
## [1] 11500 26397 17413
## [1] 11600 26790 17413
## [1] 11700 27106 17413
## [1] 11800 27494 17413
## [1] 11900 27959 17413
## [1] 12000 28351 17413
## [1] 12100 28892 17413
## [1] 12200 29303 17413## [1] 12300 29551 17413
## [1] 12400 29721 17413## [1] 12500 29857 17413
## [1] 12600 30024 17413
## [1] 12700 30174 17413
## [1] 12800 30333 17413
## [1] 12900 30505 17413
## [1] 13000 30643 17413
## [1] 13100 30827 17413
## [1] 13200 30971 17413
## [1] 13300 31118 17413
## [1] 13400 31311 17413
## [1] 13500 31458 17413
## [1] 13600 31658 17413
## [1] 13700 31811 17413
## [1] 13800 32021 17413
## [1] 13900 32177 17413
## [1] 14000 32392 17413
## [1] 14100 32552 17413
## [1] 14200 32775 17413
## [1] 14300 32941 17413
## [1] 14400 33112 17413
## [1] 14500 33355 17413
## [1] 14600 33534 17413
## [1] 14700 33790 17413
## [1] 14800 33968 17413
## [1] 14900 34146 17413
## [1] 15000 34400 17413
## [1] 15100 34577 17413
## [1] 15200 34753 17413
## [1] 15300 35004 17413
## [1] 15400 35179 17413
## [1] 15500 35353 17413
## [1] 15600 35602 17413
## [1] 15700 35777 17413
## [1] 15800 35952 17413
## [1] 15900 36129 17413
## [1] 16000 36377 17413
## [1] 16100 36553 17413
## [1] 16200 36728 17413
## [1] 16300 36975 17413
## [1] 16400 37148 17413
## [1] 16500 37391 17413
## [1] 16600 37635 17413
## [1] 16700 37877 17413
## [1] 16800 38188 17413
## [1] 16900 38498 17413
## [1] 17000 38808 17413
## [1] 17100 39119 17413
## [1] 17200 39428 17413
## [1] 17300 39759 17413
## [1] 17400 40107 17413NN.sigma_today_1800_S <- loop_sigma_D(A[S_hem,], B[S_hem,], C, "_today_1800")## [1] 100 140 22249
## [1] 200 298 22249
## [1] 300 459 22249
## [1] 400 622 22249
## [1] 500 765 22249
## [1] 600 909 22249
## [1] 700 1052 22249
## [1] 800 1195 22249
## [1] 900 1345 22249
## [1] 1000 1468 22249
## [1] 1100 1619 22249
## [1] 1200 1763 22249
## [1] 1300 1923 22249
## [1] 1400 2093 22249
## [1] 1500 2228 22249
## [1] 1600 2401 22249
## [1] 1700 2539 22249
## [1] 1800 2718 22249
## [1] 1900 2858 22249
## [1] 2000 3040 22249
## [1] 2100 3182 22249
## [1] 2200 3372 22249
## [1] 2300 3518 22249
## [1] 2400 3710 22249
## [1] 2500 3856 22249
## [1] 2600 4047 22249
## [1] 2700 4190 22249
## [1] 2800 4368 22249
## [1] 2900 4506 22249
## [1] 3000 4678 22249
## [1] 3100 4810 22249
## [1] 3200 4967 22249
## [1] 3300 5092 22249
## [1] 3400 5240 22249
## [1] 3500 5367 22249
## [1] 3600 5531 22249
## [1] 3700 5666 22249
## [1] 3800 5834 22249
## [1] 3900 5969 22249
## [1] 4000 6141 22249
## [1] 4100 6278 22249
## [1] 4200 6451 22249
## [1] 4300 6586 22249
## [1] 4400 6753 22249
## [1] 4500 6882 22249
## [1] 4600 7031 22249
## [1] 4700 7149 22249
## [1] 4800 7275 22249
## [1] 4900 7388 22249
## [1] 5000 7514 22249
## [1] 5100 7640 22249
## [1] 5200 7753 22249
## [1] 5300 7882 22249
## [1] 5400 8010 22249
## [1] 5500 8139 22249## [1] 5600 8256 22249
## [1] 5700 8404 22249
## [1] 5800 8538 22249
## [1] 5900 8680 22249
## [1] 6000 8864 22249
## [1] 6100 9073 22249
## [1] 6200 9259 22249
## [1] 6300 9501 22249
## [1] 6400 9758 22249
## [1] 6500 10019 22249
## [1] 6600 10336 22249
## [1] 6700 10614 22249
## [1] 6800 10899 22249
## [1] 6900 11128 22249
## [1] 7000 11363 22249
## [1] 7100 11532 22249
## [1] 7200 11769 22249
## [1] 7300 11938 22249
## [1] 7400 12178 22249
## [1] 7500 12350 22249
## [1] 7600 12523 22249
## [1] 7700 12772 22249
## [1] 7800 12948 22249
## [1] 7900 13203 22249
## [1] 8000 13382 22249
## [1] 8100 13561 22249
## [1] 8200 13824 22249
## [1] 8300 14012 22249
## [1] 8400 14279 22249
## [1] 8500 14462 22249
## [1] 8600 14649 22249
## [1] 8700 14915 22249
## [1] 8800 15102 22249
## [1] 8900 15289 22249
## [1] 9000 15556 22249
## [1] 9100 15740 22249
## [1] 9200 15926 22249
## [1] 9300 16112 22249
## [1] 9400 16386 22249
## [1] 9500 16573 22249
## [1] 9600 16760 22249
## [1] 9700 16947 22249
## [1] 9800 17224 22249
## [1] 9900 17414 22249
## [1] 10000 17604 22249
## [1] 10100 17881 22249
## [1] 10200 18065 22249
## [1] 10300 18247 22249
## [1] 10400 18428 22249
## [1] 10500 18688 22249
## [1] 10600 18867 22249
## [1] 10700 19044 22249
## [1] 10800 19298 22249
## [1] 10900 19475 22249
## [1] 11000 19652 22249
## [1] 11100 19828 22249
## [1] 11200 20081 22249
## [1] 11300 20258 22249
## [1] 11400 20436 22249
## [1] 11500 20694 22249
## [1] 11600 20874 22249
## [1] 11700 21054 22249
## [1] 11800 21314 22249
## [1] 11900 21494 22249
## [1] 12000 21674 22249
## [1] 12100 21934 22249
## [1] 12200 22115 22249
## [1] 12300 22295 22249
## [1] 12400 22553 22249
## [1] 12500 22731 22249
## [1] 12600 22909 22249
## [1] 12700 23161 22249
## [1] 12800 23335 22249
## [1] 12900 23508 22249
## [1] 13000 23753 22249
## [1] 13100 23923 22249
## [1] 13200 24092 22249
## [1] 13300 24326 22249
## [1] 13400 24485 22249
## [1] 13500 24690 22249
## [1] 13600 24839 22249
## [1] 13700 24987 22249
## [1] 13800 25178 22249
## [1] 13900 25322 22249
## [1] 14000 25465 22249
## [1] 14100 25646 22249
## [1] 14200 25785 22249
## [1] 14300 25923 22249
## [1] 14400 26095 22249
## [1] 14500 26228 22249
## [1] 14600 26360 22249
## [1] 14700 26522 22249
## [1] 14800 26652 22249
## [1] 14900 26811 22249
## [1] 15000 26938 22249
## [1] 15100 27065 22249
## [1] 15200 27218 22249
## [1] 15300 27344 22249
## [1] 15400 27470 22249
## [1] 15500 27619 22249
## [1] 15600 27743 22249
## [1] 15700 27890 22249
## [1] 15800 28011 22249
## [1] 15900 28132 22249
## [1] 16000 28274 22249
## [1] 16100 28394 22249
## [1] 16200 28514 22249
## [1] 16300 28648 22249
## [1] 16400 28764 22249
## [1] 16500 28898 22249
## [1] 16600 29022 22249
## [1] 16700 29172 22249
## [1] 16800 29327 22249
## [1] 16900 29498 22249
## [1] 17000 29708 22249
## [1] 17100 30060 22249
## [1] 17200 30380 22249
## [1] 17300 30680 22249
## [1] 17400 30923 22249
## [1] 17500 31234 22249
## [1] 17600 31474 22249
## [1] 17700 31716 22249
## [1] 17800 31954 22249
## [1] 17900 32192 22249
## [1] 18000 32359 22249
## [1] 18100 32598 22249
## [1] 18200 32844 22249
## [1] 18300 33019 22249
## [1] 18400 33266 22249
## [1] 18500 33439 22249
## [1] 18600 33616 22249
## [1] 18700 33863 22249
## [1] 18800 34038 22249
## [1] 18900 34213 22249
## [1] 19000 34463 22249
## [1] 19100 34638 22249
## [1] 19200 34814 22249
## [1] 19300 35066 22249
## [1] 19400 35242 22249
## [1] 19500 35419 22249
## [1] 19600 35678 22249
## [1] 19700 35859 22249
## [1] 19800 36042 22249
## [1] 19900 36312 22249
## [1] 20000 36490 22249
## [1] 20100 36664 22249
## [1] 20200 36908 22249
## [1] 20300 37078 22249
## [1] 20400 37310 22249
## [1] 20500 37464 22249
## [1] 20600 37612 22249
## [1] 20700 37801 22249
## [1] 20800 37943 22249
## [1] 20900 38122 22249
## [1] 21000 38258 22249
## [1] 21100 38428 22249
## [1] 21200 38563 22249
## [1] 21300 38697 22249
## [1] 21400 38864 22249
## [1] 21500 38997 22249
## [1] 21600 39161 22249
## [1] 21700 39292 22249
## [1] 21800 39446 22249
## [1] 21900 39572 22249
## [1] 22000 39721 22249
## [1] 22100 39866 22249
## [1] 22200 40030 22249NN.sigma_today_1800 <- rbind(NN.sigma_today_1800_N,
NN.sigma_today_1800_S)
head(NN.sigma_today_1800)## No NN.sigma_today_1800 NN.station_today_1800 NN.Mdist_today_1800
## 1 36 0.836709022 31931 1.3486467
## 2 37 0.894133247 794 1.4077485
## 3 38 0.661532267 861 1.1633920
## 4 39 0.265489647 39602 0.6854521
## 5 40 0.006257156 39603 0.1000496
## 6 41 0.036266868 39426 0.2423104
## numPCs_today_1800
## 1 2
## 2 2
## 3 2
## 4 2
## 5 2
## 6 2final_dat2 <- full_join(final_dat, NN.sigma_today_1800)## Joining, by = "No"head(final_dat2)## No lat long NN.sigma_1800_today NN.station_1800_today NN.Mdist_1800_today
## 1 1 -69.5 20.5 0.005340623 40066 0.0924152
## 2 2 -68.5 20.5 0.007501424 40009 0.1095737
## 3 3 -67.5 20.5 0.009629160 40067 0.1241975
## 4 4 -66.5 20.5 0.015490226 40009 0.1577081
## 5 5 -65.5 20.5 0.024201364 40066 0.1974682
## 6 6 -64.5 20.5 0.049860481 40066 0.2848815
## numPCs_1800_today NN.sigma_today_1800 NN.station_today_1800
## 1 2 0.019575390 26505
## 2 2 0.020981196 29981
## 3 2 0.014089519 26082
## 4 2 0.006465863 26082
## 5 2 0.011645698 25749
## 6 2 0.005112661 25973
## NN.Mdist_today_1800 numPCs_today_1800
## 1 0.17743251 2
## 2 0.18374455 2
## 3 0.15036689 2
## 4 0.10170870 2
## 5 0.13663930 2
## 6 0.09041724 2dim(final_dat2)## [1] 39662 11Plot_nonInt(final_dat2$lat, final_dat2$long,
final_dat2$NN.sigma_today_1800, world, "sigD dis 2000")
ggplot(final_dat2) + geom_point(aes(x=NN.sigma_1800_today,
y=NN.sigma_today_1800,
color=lat), alpha=0.5) +
scale_color_gradient2(low="red", high="blue", mid="grey") +
theme_classic() + xlab("Novel climates today from 1800") +
ylab("Disappearing 1800 climates")
# The NN.sigma_1800_disappear represents the novelty of 1800
# compared to today (e.g disappearing climates)
# NN.dist_1800_disappear represents the Mahalanobis distance from the climate
# in 1800 to it's nearest neighbor in all of today's data
# No--> NN.station_1800_disappear represents the the station that was queried 1800 (B)
# and it's nearest data in the global reference NN.station_1800_disappear (2000, A)Add the info for lat/long of nearest neighbors (1800 to today)
stationInfo2 <- stationInfo
names(stationInfo2) <- paste0(names(stationInfo2),"_1800_today")
head(stationInfo2)## No_1800_today lat_1800_today long_1800_today
## 1 1 -69.5 20.5
## 2 2 -68.5 20.5
## 3 3 -67.5 20.5
## 4 4 -66.5 20.5
## 5 5 -65.5 20.5
## 6 6 -64.5 20.5dim(stationInfo)## [1] 39662 3names(stationInfo2)[1] <-"NN.station_1800_today"
final_dat3 <- left_join(final_dat2, stationInfo2)## Joining, by = "NN.station_1800_today"dim(final_dat3)## [1] 39662 13head(final_dat3)## No lat long NN.sigma_1800_today NN.station_1800_today NN.Mdist_1800_today
## 1 1 -69.5 20.5 0.005340623 40066 0.0924152
## 2 2 -68.5 20.5 0.007501424 40009 0.1095737
## 3 3 -67.5 20.5 0.009629160 40067 0.1241975
## 4 4 -66.5 20.5 0.015490226 40009 0.1577081
## 5 5 -65.5 20.5 0.024201364 40066 0.1974682
## 6 6 -64.5 20.5 0.049860481 40066 0.2848815
## numPCs_1800_today NN.sigma_today_1800 NN.station_today_1800
## 1 2 0.019575390 26505
## 2 2 0.020981196 29981
## 3 2 0.014089519 26082
## 4 2 0.006465863 26082
## 5 2 0.011645698 25749
## 6 2 0.005112661 25973
## NN.Mdist_today_1800 numPCs_today_1800 lat_1800_today long_1800_today
## 1 0.17743251 2 -69.5 379.5
## 2 0.18374455 2 -69.5 378.5
## 3 0.15036689 2 -68.5 379.5
## 4 0.10170870 2 -69.5 378.5
## 5 0.13663930 2 -69.5 379.5
## 6 0.09041724 2 -69.5 379.5tail(final_dat3)## No lat long NN.sigma_1800_today NN.station_1800_today
## 39657 40115 84.5 379.5 0.0033782225 37019
## 39658 40116 85.5 379.5 0.0041649559 37769
## 39659 40117 86.5 379.5 0.0001037297 3229
## 39660 40118 87.5 379.5 0.0047175381 7551
## 39661 40119 88.5 379.5 0.0377153031 7551
## 39662 40120 89.5 379.5 0.0646634655 7551
## NN.Mdist_1800_today numPCs_1800_today NN.sigma_today_1800
## 39657 0.07347198 2 1.728258
## 39658 0.08159263 2 1.931682
## 39659 0.01286606 2 2.170640
## 39660 0.08684629 2 2.350300
## 39661 0.24717272 2 2.424026
## 39662 0.32537576 2 2.455741
## NN.station_today_1800 NN.Mdist_today_1800 numPCs_today_1800
## 39657 613 2.226041 2
## 39658 613 2.420731 2
## 39659 613 2.648753 2
## 39660 554 2.819972 2
## 39661 3223 2.890216 2
## 39662 3336 2.920432 2
## lat_1800_today long_1800_today
## 39657 86.5 347.5
## 39658 86.5 353.5
## 39659 88.5 59.5
## 39660 88.5 103.5
## 39661 88.5 103.5
## 39662 88.5 103.5# Should both be 0
sum(!complete.cases(final_dat2))## [1] 0sum(!complete.cases(final_dat3))## [1] 0head(final_dat3[!complete.cases(final_dat3),])## [1] No lat long
## [4] NN.sigma_1800_today NN.station_1800_today NN.Mdist_1800_today
## [7] numPCs_1800_today NN.sigma_today_1800 NN.station_today_1800
## [10] NN.Mdist_today_1800 numPCs_today_1800 lat_1800_today
## [13] long_1800_today
## <0 rows> (or 0-length row.names)summary(final_dat3, NN.station_today_1800)## No lat long NN.sigma_1800_today
## Min. : 1 Min. :-78.500 Min. : 20.5 Min. :0.0000002
## 1st Qu.:10100 1st Qu.:-43.500 1st Qu.:138.5 1st Qu.:0.0037303
## Median :20110 Median : -9.500 Median :207.5 Median :0.0778035
## Mean :20095 Mean : -1.731 Mean :208.8 Mean :0.4287141
## 3rd Qu.:30087 3rd Qu.: 34.500 3rd Qu.:295.5 3rd Qu.:0.6417823
## Max. :40120 Max. : 89.500 Max. :379.5 Max. :3.1374782
## NN.station_1800_today NN.Mdist_1800_today numPCs_1800_today
## Min. : 44 Min. :0.000572 Min. :2.000
## 1st Qu.:20665 1st Qu.:0.077212 1st Qu.:2.000
## Median :25736 Median :0.357832 Median :2.000
## Mean :24568 Mean :0.705606 Mean :2.002
## 3rd Qu.:29342 3rd Qu.:1.142396 3rd Qu.:2.000
## Max. :40074 Max. :3.611830 Max. :3.000
## NN.sigma_today_1800 NN.station_today_1800 NN.Mdist_today_1800
## Min. :0.00000 Min. : 39 Min. :0.000314
## 1st Qu.:0.03089 1st Qu.:17767 1st Qu.:0.224749
## Median :0.83348 Median :30424 Median :1.345494
## Mean :0.92699 Mean :24144 Mean :1.296300
## 3rd Qu.:1.60850 3rd Qu.:31698 3rd Qu.:2.112351
## Max. :3.44851 Max. :40119 Max. :3.868075
## numPCs_today_1800 lat_1800_today long_1800_today
## Min. :2.000 Min. :-77.500 Min. : 20.5
## 1st Qu.:2.000 1st Qu.:-43.500 1st Qu.:210.5
## Median :2.000 Median : -2.500 Median :246.5
## Mean :2.002 Mean : -2.377 Mean :244.0
## 3rd Qu.:2.000 3rd Qu.: 25.500 3rd Qu.:284.5
## Max. :3.000 Max. : 89.500 Max. :379.5final_dat3$latshift_1800_today <- NA
condN <- which(final_dat3$lat>0)
final_dat3$latshift_1800_today[condN] <- final_dat3$lat[condN] -
final_dat3$lat_1800_today[condN]
condS <- which(final_dat3$lat<0)
final_dat3$latshift[condS] <- -1*(final_dat3$lat[condS] -
final_dat3$lat_1800_today[condS])
hist(final_dat3$latshift_1800_today)
#stationInfo3 <- stationInfo
#names(stationInfo3) <- paste0(names(stationInfo3),"_today_1800")
#names(stationInfo3)[1] <- "NN.station_today_1800"
#final_dat4 <- left_join(final_dat3, stationInfo3)
#dim(final_dat4)
#head(final_dat4)
# cond <- which(final_dat4$NN.sigma_1800_today>5 | final_dat4$NN.sigma_1800_disappear>5)
# length(cond)
# final_dat4[cond,]ggplot(final_dat3) + geom_point(aes(y=lat,
x=lat_1800_today,
color=NN.sigma_1800_today), alpha=0.2) +
scale_color_gradient2(low="red", high="blue", mid="grey") +
theme_classic() + ylab("Latitude in 2000") +
xlab("Latitude of nearest neighbor in 1800") + geom_abline(intercept=0,slope=1) 
# This plot makes more sense in the way we typically think about it.
# If we consider a location in 2000, where was the climate it came from in 1800?
write.csv(final_dat3, paste0(results_dir, "1800base_2000_novelty_NShemrestrict.csv"))
final_dat4 <- final_dat3# ggplot(final_dat4) + geom_point(aes(x=lat_1800_disappear,
# y=lat,
# color=NN.sigma_1800_disappear), alpha=0.2) +
# scale_color_gradient2(low="red", high="blue", mid="grey") +
# theme_classic() + ylab("Latitude in 1800") +
# xlab("Latitude of nearest neighbor in 2000") + geom_abline(intercept=0,slope=1) +
# geom_abline(intercept=0,slope=-1)
# this is weird
#n <- final_dat4 %>% filter(lat_1800_disappear > 0 & lat_1800_disappear < 10)
#Plot_nonInt(n$lat_1800_disappear, n$long_1800_disappear,
# n$NN.sigma_1800_disappear, world, "sigma dis.")
#hist(final_dat4$lat_1800_disappear, breaks=seq(-90,90, by=1))
#hist(final_dat4$lat_1800_today, breaks=seq(-90,90, by=1))
#hist(final_dat4$lat, breaks=seq(-90,90, by=1))Today analog to 2100 RCP 8.5 (Novelty)
identical(norm_2000$No, norm_2100_8.5$No)## [1] TRUEA1 <- norm_2000
# 1970-2000 climate normals
head(A1)## No SST Arag pH
## 1 1 -1.4479254 0.1799026 8.035921
## 2 2 -1.3696480 0.1741309 8.045892
## 3 3 -1.2845583 0.1683092 8.053425
## 4 4 -1.1916660 0.1624748 8.058488
## 5 5 -1.0895625 0.1567233 8.061231
## 6 6 -0.9764333 0.1512065 8.062000dim(A1)## [1] 39662 4B1 <- norm_2100_8.5
# 2070-2100
head(B1)## No SST Arag pH
## 1 1 -0.55466154 -0.09210795 7.726612
## 2 2 -0.45603873 -0.09383454 7.740267
## 3 3 -0.34939250 -0.09577437 7.751387
## 4 4 -0.23345842 -0.09793483 7.759896
## 5 5 -0.10654381 -0.10027507 7.765894
## 6 6 0.03343595 -0.10269243 7.769621dim(B1)## [1] 39662 4# same columns as A
# sanity check to make sure stations in right order
identical(A1$No, B1$No) # should be true## [1] TRUEidentical(A1$No, sort(A1$No)) # should also be true## [1] TRUE# check list 18906,18952, 29736, 29789, 8193, 29319,
# calc_sigma_D(A1, B1, C, 18906, "_today_2100_8.5")
# calc_sigma_D(A1, B1, C, 18952, "_today_2100_8.5")
# calc_sigma_D(A1, B1, C, 29736, "_today_2100_8.5")
# calc_sigma_D(A1, B1, C, 29789, "_today_2100_8.5")
# calc_sigma_D(A1, B1, C, 8193, "_today_2100_8.5")
# calc_sigma_D(A1, B1, C, 29319, "_today_2100_8.5")
# C is the same defined above
NN.sigma_today_2100_8.5_N <- loop_sigma_D(A1[N_hem,], B1[N_hem,], C, "_today_2100_8.5")## [1] 100 317 17413
## [1] 200 610 17413
## [1] 300 934 17413
## [1] 400 1247 17413
## [1] 500 1683 17413
## [1] 600 2053 17413
## [1] 700 2356 17413
## [1] 800 2660 17413
## [1] 900 2971 17413
## [1] 1000 3213 17413
## [1] 1100 3452 17413
## [1] 1200 3759 17413
## [1] 1300 3996 17413
## [1] 1400 4236 17413
## [1] 1500 4540 17413
## [1] 1600 4847 17413
## [1] 1700 5220 17413
## [1] 1800 5530 17413
## [1] 1900 5830 17413
## [1] 2000 6129 17413
## [1] 2100 6427 17413
## [1] 2200 6726 17413
## [1] 2300 7099 17413
## [1] 2400 7694 17413
## [1] 2500 8267 17413
## [1] 2600 8629 17413
## [1] 2700 8887 17413
## [1] 2800 9069 17413
## [1] 2900 9286 17413
## [1] 3000 9454 17413
## [1] 3100 9622 17413
## [1] 3200 9790 17413
## [1] 3300 9954 17413
## [1] 3400 10119 17413
## [1] 3500 10287 17413
## [1] 3600 10458 17413
## [1] 3700 10592 17413
## [1] 3800 10764 17413
## [1] 3900 10945 17413
## [1] 4000 11095 17413
## [1] 4100 11309 17413
## [1] 4200 11469 17413
## [1] 4300 11702 17413
## [1] 4400 11871 17413
## [1] 4500 12109 17413
## [1] 4600 12278 17413
## [1] 4700 12450 17413
## [1] 4800 12692 17413
## [1] 4900 12863 17413
## [1] 5000 13035 17413
## [1] 5100 13280 17413
## [1] 5200 13453 17413
## [1] 5300 13627 17413
## [1] 5400 13878 17413
## [1] 5500 14056 17413
## [1] 5600 14235 17413
## [1] 5700 14414 17413
## [1] 5800 14670 17413
## [1] 5900 14848 17413
## [1] 6000 15026 17413
## [1] 6100 15204 17413
## [1] 6200 15384 17413
## [1] 6300 15641 17413
## [1] 6400 15819 17413
## [1] 6500 15998 17413
## [1] 6600 16176 17413
## [1] 6700 16356 17413
## [1] 6800 16534 17413
## [1] 6900 16791 17413
## [1] 7000 16970 17413
## [1] 7100 17149 17413
## [1] 7200 17328 17413
## [1] 7300 17506 17413
## [1] 7400 17684 17413
## [1] 7500 17862 17413
## [1] 7600 18119 17413
## [1] 7700 18297 17413
## [1] 7800 18476 17413
## [1] 7900 18655 17413## [1] 8000 18912 17413## [1] 8100 19090 17413
## [1] 8200 19269 17413
## [1] 8300 19526 17413
## [1] 8400 19704 17413
## [1] 8500 19881 17413
## [1] 8600 20058 17413
## [1] 8700 20312 17413
## [1] 8800 20488 17413
## [1] 8900 20664 17413
## [1] 9000 20916 17413
## [1] 9100 21091 17413
## [1] 9200 21266 17413
## [1] 9300 21441 17413
## [1] 9400 21691 17413
## [1] 9500 21866 17413
## [1] 9600 22041 17413
## [1] 9700 22218 17413
## [1] 9800 22468 17413
## [1] 9900 22643 17413
## [1] 10000 22817 17413
## [1] 10100 22992 17413
## [1] 10200 23240 17413
## [1] 10300 23414 17413
## [1] 10400 23662 17413
## [1] 10500 23836 17413
## [1] 10600 24009 17413
## [1] 10700 24256 17413
## [1] 10800 24430 17413
## [1] 10900 24679 17413
## [1] 11000 24927 17413
## [1] 11100 25250 17413
## [1] 11200 25500 17413
## [1] 11300 25823 17413
## [1] 11400 26075 17413
## [1] 11500 26397 17413
## [1] 11600 26790 17413
## [1] 11700 27106 17413
## [1] 11800 27494 17413
## [1] 11900 27959 17413
## [1] 12000 28351 17413
## [1] 12100 28892 17413
## [1] 12200 29303 17413## [1] 12300 29551 17413
## [1] 12400 29721 17413## [1] 12500 29857 17413
## [1] 12600 30024 17413
## [1] 12700 30174 17413
## [1] 12800 30333 17413
## [1] 12900 30505 17413
## [1] 13000 30643 17413
## [1] 13100 30827 17413
## [1] 13200 30971 17413
## [1] 13300 31118 17413
## [1] 13400 31311 17413
## [1] 13500 31458 17413
## [1] 13600 31658 17413
## [1] 13700 31811 17413
## [1] 13800 32021 17413
## [1] 13900 32177 17413
## [1] 14000 32392 17413
## [1] 14100 32552 17413
## [1] 14200 32775 17413
## [1] 14300 32941 17413
## [1] 14400 33112 17413
## [1] 14500 33355 17413
## [1] 14600 33534 17413
## [1] 14700 33790 17413
## [1] 14800 33968 17413
## [1] 14900 34146 17413
## [1] 15000 34400 17413
## [1] 15100 34577 17413
## [1] 15200 34753 17413
## [1] 15300 35004 17413
## [1] 15400 35179 17413
## [1] 15500 35353 17413
## [1] 15600 35602 17413
## [1] 15700 35777 17413
## [1] 15800 35952 17413
## [1] 15900 36129 17413
## [1] 16000 36377 17413
## [1] 16100 36553 17413
## [1] 16200 36728 17413
## [1] 16300 36975 17413
## [1] 16400 37148 17413
## [1] 16500 37391 17413
## [1] 16600 37635 17413
## [1] 16700 37877 17413
## [1] 16800 38188 17413
## [1] 16900 38498 17413
## [1] 17000 38808 17413
## [1] 17100 39119 17413
## [1] 17200 39428 17413
## [1] 17300 39759 17413
## [1] 17400 40107 17413# C is the same defined above
NN.sigma_today_2100_8.5_S <- loop_sigma_D(A1[S_hem,], B1[S_hem,], C, "_today_2100_8.5")## [1] 100 140 22249
## [1] 200 298 22249
## [1] 300 459 22249
## [1] 400 622 22249
## [1] 500 765 22249
## [1] 600 909 22249
## [1] 700 1052 22249
## [1] 800 1195 22249
## [1] 900 1345 22249
## [1] 1000 1468 22249
## [1] 1100 1619 22249
## [1] 1200 1763 22249
## [1] 1300 1923 22249
## [1] 1400 2093 22249
## [1] 1500 2228 22249
## [1] 1600 2401 22249
## [1] 1700 2539 22249
## [1] 1800 2718 22249
## [1] 1900 2858 22249
## [1] 2000 3040 22249
## [1] 2100 3182 22249
## [1] 2200 3372 22249
## [1] 2300 3518 22249
## [1] 2400 3710 22249
## [1] 2500 3856 22249
## [1] 2600 4047 22249
## [1] 2700 4190 22249
## [1] 2800 4368 22249
## [1] 2900 4506 22249
## [1] 3000 4678 22249
## [1] 3100 4810 22249
## [1] 3200 4967 22249
## [1] 3300 5092 22249
## [1] 3400 5240 22249
## [1] 3500 5367 22249
## [1] 3600 5531 22249
## [1] 3700 5666 22249
## [1] 3800 5834 22249
## [1] 3900 5969 22249
## [1] 4000 6141 22249
## [1] 4100 6278 22249
## [1] 4200 6451 22249
## [1] 4300 6586 22249
## [1] 4400 6753 22249
## [1] 4500 6882 22249
## [1] 4600 7031 22249
## [1] 4700 7149 22249
## [1] 4800 7275 22249
## [1] 4900 7388 22249
## [1] 5000 7514 22249
## [1] 5100 7640 22249
## [1] 5200 7753 22249
## [1] 5300 7882 22249
## [1] 5400 8010 22249
## [1] 5500 8139 22249## [1] 5600 8256 22249
## [1] 5700 8404 22249
## [1] 5800 8538 22249
## [1] 5900 8680 22249
## [1] 6000 8864 22249
## [1] 6100 9073 22249
## [1] 6200 9259 22249
## [1] 6300 9501 22249
## [1] 6400 9758 22249
## [1] 6500 10019 22249
## [1] 6600 10336 22249
## [1] 6700 10614 22249
## [1] 6800 10899 22249
## [1] 6900 11128 22249
## [1] 7000 11363 22249
## [1] 7100 11532 22249
## [1] 7200 11769 22249
## [1] 7300 11938 22249
## [1] 7400 12178 22249
## [1] 7500 12350 22249
## [1] 7600 12523 22249
## [1] 7700 12772 22249
## [1] 7800 12948 22249
## [1] 7900 13203 22249
## [1] 8000 13382 22249
## [1] 8100 13561 22249
## [1] 8200 13824 22249
## [1] 8300 14012 22249
## [1] 8400 14279 22249
## [1] 8500 14462 22249
## [1] 8600 14649 22249
## [1] 8700 14915 22249
## [1] 8800 15102 22249
## [1] 8900 15289 22249
## [1] 9000 15556 22249
## [1] 9100 15740 22249
## [1] 9200 15926 22249
## [1] 9300 16112 22249
## [1] 9400 16386 22249
## [1] 9500 16573 22249
## [1] 9600 16760 22249
## [1] 9700 16947 22249
## [1] 9800 17224 22249
## [1] 9900 17414 22249
## [1] 10000 17604 22249
## [1] 10100 17881 22249
## [1] 10200 18065 22249
## [1] 10300 18247 22249
## [1] 10400 18428 22249
## [1] 10500 18688 22249
## [1] 10600 18867 22249
## [1] 10700 19044 22249
## [1] 10800 19298 22249
## [1] 10900 19475 22249
## [1] 11000 19652 22249
## [1] 11100 19828 22249
## [1] 11200 20081 22249
## [1] 11300 20258 22249
## [1] 11400 20436 22249
## [1] 11500 20694 22249
## [1] 11600 20874 22249
## [1] 11700 21054 22249
## [1] 11800 21314 22249
## [1] 11900 21494 22249
## [1] 12000 21674 22249
## [1] 12100 21934 22249
## [1] 12200 22115 22249
## [1] 12300 22295 22249
## [1] 12400 22553 22249
## [1] 12500 22731 22249
## [1] 12600 22909 22249
## [1] 12700 23161 22249
## [1] 12800 23335 22249
## [1] 12900 23508 22249
## [1] 13000 23753 22249
## [1] 13100 23923 22249
## [1] 13200 24092 22249
## [1] 13300 24326 22249
## [1] 13400 24485 22249
## [1] 13500 24690 22249
## [1] 13600 24839 22249
## [1] 13700 24987 22249
## [1] 13800 25178 22249
## [1] 13900 25322 22249
## [1] 14000 25465 22249
## [1] 14100 25646 22249
## [1] 14200 25785 22249
## [1] 14300 25923 22249
## [1] 14400 26095 22249
## [1] 14500 26228 22249
## [1] 14600 26360 22249
## [1] 14700 26522 22249
## [1] 14800 26652 22249
## [1] 14900 26811 22249
## [1] 15000 26938 22249
## [1] 15100 27065 22249
## [1] 15200 27218 22249
## [1] 15300 27344 22249
## [1] 15400 27470 22249
## [1] 15500 27619 22249
## [1] 15600 27743 22249
## [1] 15700 27890 22249
## [1] 15800 28011 22249
## [1] 15900 28132 22249
## [1] 16000 28274 22249
## [1] 16100 28394 22249
## [1] 16200 28514 22249
## [1] 16300 28648 22249
## [1] 16400 28764 22249
## [1] 16500 28898 22249
## [1] 16600 29022 22249
## [1] 16700 29172 22249
## [1] 16800 29327 22249
## [1] 16900 29498 22249
## [1] 17000 29708 22249
## [1] 17100 30060 22249
## [1] 17200 30380 22249
## [1] 17300 30680 22249
## [1] 17400 30923 22249
## [1] 17500 31234 22249
## [1] 17600 31474 22249
## [1] 17700 31716 22249
## [1] 17800 31954 22249
## [1] 17900 32192 22249
## [1] 18000 32359 22249
## [1] 18100 32598 22249
## [1] 18200 32844 22249
## [1] 18300 33019 22249
## [1] 18400 33266 22249
## [1] 18500 33439 22249
## [1] 18600 33616 22249
## [1] 18700 33863 22249
## [1] 18800 34038 22249
## [1] 18900 34213 22249
## [1] 19000 34463 22249
## [1] 19100 34638 22249
## [1] 19200 34814 22249
## [1] 19300 35066 22249
## [1] 19400 35242 22249
## [1] 19500 35419 22249
## [1] 19600 35678 22249
## [1] 19700 35859 22249
## [1] 19800 36042 22249
## [1] 19900 36312 22249
## [1] 20000 36490 22249
## [1] 20100 36664 22249
## [1] 20200 36908 22249
## [1] 20300 37078 22249
## [1] 20400 37310 22249
## [1] 20500 37464 22249
## [1] 20600 37612 22249
## [1] 20700 37801 22249
## [1] 20800 37943 22249
## [1] 20900 38122 22249
## [1] 21000 38258 22249
## [1] 21100 38428 22249
## [1] 21200 38563 22249
## [1] 21300 38697 22249
## [1] 21400 38864 22249
## [1] 21500 38997 22249
## [1] 21600 39161 22249
## [1] 21700 39292 22249
## [1] 21800 39446 22249
## [1] 21900 39572 22249
## [1] 22000 39721 22249
## [1] 22100 39866 22249
## [1] 22200 40030 22249NN.sigma_today_2100_8.5 <- rbind(NN.sigma_today_2100_8.5_N,
NN.sigma_today_2100_8.5_S)
which(is.infinite(NN.sigma_today_2100_8.5$NN.sigma_today_2100_8.5))## integer(0)which(is.na(NN.sigma_today_2100_8.5$NN.sigma_today_2100_8.5))## integer(0)final_dat5 <- full_join(NN.sigma_today_2100_8.5, final_dat4, by="No")
head(final_dat5)## No NN.sigma_today_2100_8.5 NN.station_today_2100_8.5 NN.Mdist_today_2100_8.5
## 1 36 0.3501681236 29638 0.79989070
## 2 37 0.1043489359 29695 0.41656957
## 3 38 0.0015820951 29638 0.05026184
## 4 39 0.0006816715 28978 0.03298615
## 5 40 0.1188184652 28440 0.44577407
## 6 41 0.5463419465 28254 1.03579199
## numPCs_today_2100_8.5 lat long NN.sigma_1800_today NN.station_1800_today
## 1 2 70.5 20.5 0.033725924 30751
## 2 2 71.5 20.5 0.015960753 30751
## 3 2 72.5 20.5 0.010016616 31111
## 4 2 73.5 20.5 0.001600598 14563
## 5 2 74.5 20.5 0.018324556 1958
## 6 2 75.5 20.5 0.004808260 17177
## NN.Mdist_1800_today numPCs_1800_today NN.sigma_today_1800
## 1 0.23355016 2 0.836709022
## 2 0.16010043 2 0.894133247
## 3 0.12668130 2 0.661532267
## 4 0.05055509 2 0.265489647
## 5 0.17162744 2 0.006257156
## 6 0.08767895 2 0.036266868
## NN.station_today_1800 NN.Mdist_today_1800 numPCs_today_1800 lat_1800_today
## 1 31931 1.3486467 2 44.5
## 2 794 1.4077485 2 44.5
## 3 861 1.1633920 2 82.5
## 4 39602 0.6854521 2 53.5
## 5 39603 0.1000496 2 69.5
## 6 39426 0.2423104 2 61.5
## long_1800_today latshift_1800_today latshift
## 1 302.5 26 26
## 2 302.5 27 27
## 3 305.5 -10 -10
## 4 172.5 20 20
## 5 47.5 5 5
## 6 188.5 14 14dim(final_dat5)## [1] 39662 19stationInfo5 <- stationInfo
names(stationInfo5) <- paste0(names(stationInfo),"_today_2100_8.5")
names(stationInfo5)[1] <- "NN.station_today_2100_8.5"
head(stationInfo5)## NN.station_today_2100_8.5 lat_today_2100_8.5 long_today_2100_8.5
## 1 1 -69.5 20.5
## 2 2 -68.5 20.5
## 3 3 -67.5 20.5
## 4 4 -66.5 20.5
## 5 5 -65.5 20.5
## 6 6 -64.5 20.5final_dat6 <- left_join(final_dat5, stationInfo5)## Joining, by = "NN.station_today_2100_8.5"dim(final_dat6)## [1] 39662 21head(final_dat6)## No NN.sigma_today_2100_8.5 NN.station_today_2100_8.5 NN.Mdist_today_2100_8.5
## 1 36 0.3501681236 29638 0.79989070
## 2 37 0.1043489359 29695 0.41656957
## 3 38 0.0015820951 29638 0.05026184
## 4 39 0.0006816715 28978 0.03298615
## 5 40 0.1188184652 28440 0.44577407
## 6 41 0.5463419465 28254 1.03579199
## numPCs_today_2100_8.5 lat long NN.sigma_1800_today NN.station_1800_today
## 1 2 70.5 20.5 0.033725924 30751
## 2 2 71.5 20.5 0.015960753 30751
## 3 2 72.5 20.5 0.010016616 31111
## 4 2 73.5 20.5 0.001600598 14563
## 5 2 74.5 20.5 0.018324556 1958
## 6 2 75.5 20.5 0.004808260 17177
## NN.Mdist_1800_today numPCs_1800_today NN.sigma_today_1800
## 1 0.23355016 2 0.836709022
## 2 0.16010043 2 0.894133247
## 3 0.12668130 2 0.661532267
## 4 0.05055509 2 0.265489647
## 5 0.17162744 2 0.006257156
## 6 0.08767895 2 0.036266868
## NN.station_today_1800 NN.Mdist_today_1800 numPCs_today_1800 lat_1800_today
## 1 31931 1.3486467 2 44.5
## 2 794 1.4077485 2 44.5
## 3 861 1.1633920 2 82.5
## 4 39602 0.6854521 2 53.5
## 5 39603 0.1000496 2 69.5
## 6 39426 0.2423104 2 61.5
## long_1800_today latshift_1800_today latshift lat_today_2100_8.5
## 1 302.5 26 26 89.5
## 2 302.5 27 27 89.5
## 3 305.5 -10 -10 89.5
## 4 172.5 20 20 85.5
## 5 47.5 5 5 82.5
## 6 188.5 14 14 82.5
## long_today_2100_8.5
## 1 287.5
## 2 288.5
## 3 287.5
## 4 279.5
## 5 273.5
## 6 271.5# Visualize
Plot_nonInt(final_dat6$lat, final_dat6$long,
final_dat6$NN.sigma_today_2100_8.5, world, "sigma nov 2100 8.5")
2100 RCP 8.5 to today, what are today’s climates that will disappear in 2100 RCP 8.5?
#--------------------------------
### 2100 RCP 8.5 to today, what are today's climates that ####
### will disappear in 2100 RCP 8.5?
#--------------------------------
# check list 18906,18952, 29736, 29789, 8193, 29319,
# calc_sigma_D( B1, A1, C, 18906, "_2100_8.5_today")
# calc_sigma_D( B1, A1, C, 18952, "_2100_8.5_today")
# calc_sigma_D( B1, A1, C, 29736, "_2100_8.5_today")
# calc_sigma_D( B1, A1, C, 29789, "_2100_8.5_today")
# calc_sigma_D( B1, A1, C, 8193, "_2100_8.5_today")
# calc_sigma_D( B1, A1, C, 29319, "_2100_8.5_today")
NN.sigma_2100_8.5_today_N <- loop_sigma_D( B1[N_hem,], A1[N_hem,], C, "_2100_8.5_today")## [1] 100 317 17413
## [1] 200 610 17413
## [1] 300 934 17413
## [1] 400 1247 17413
## [1] 500 1683 17413
## [1] 600 2053 17413
## [1] 700 2356 17413
## [1] 800 2660 17413
## [1] 900 2971 17413
## [1] 1000 3213 17413
## [1] 1100 3452 17413
## [1] 1200 3759 17413
## [1] 1300 3996 17413
## [1] 1400 4236 17413
## [1] 1500 4540 17413
## [1] 1600 4847 17413
## [1] 1700 5220 17413
## [1] 1800 5530 17413
## [1] 1900 5830 17413
## [1] 2000 6129 17413
## [1] 2100 6427 17413
## [1] 2200 6726 17413
## [1] 2300 7099 17413
## [1] 2400 7694 17413
## [1] 2500 8267 17413
## [1] 2600 8629 17413
## [1] 2700 8887 17413
## [1] 2800 9069 17413
## [1] 2900 9286 17413
## [1] 3000 9454 17413
## [1] 3100 9622 17413
## [1] 3200 9790 17413
## [1] 3300 9954 17413
## [1] 3400 10119 17413
## [1] 3500 10287 17413
## [1] 3600 10458 17413
## [1] 3700 10592 17413
## [1] 3800 10764 17413
## [1] 3900 10945 17413
## [1] 4000 11095 17413
## [1] 4100 11309 17413
## [1] 4200 11469 17413
## [1] 4300 11702 17413
## [1] 4400 11871 17413
## [1] 4500 12109 17413
## [1] 4600 12278 17413
## [1] 4700 12450 17413
## [1] 4800 12692 17413
## [1] 4900 12863 17413
## [1] 5000 13035 17413
## [1] 5100 13280 17413
## [1] 5200 13453 17413
## [1] 5300 13627 17413
## [1] 5400 13878 17413
## [1] 5500 14056 17413
## [1] 5600 14235 17413
## [1] 5700 14414 17413
## [1] 5800 14670 17413
## [1] 5900 14848 17413
## [1] 6000 15026 17413
## [1] 6100 15204 17413
## [1] 6200 15384 17413
## [1] 6300 15641 17413
## [1] 6400 15819 17413
## [1] 6500 15998 17413
## [1] 6600 16176 17413
## [1] 6700 16356 17413
## [1] 6800 16534 17413
## [1] 6900 16791 17413
## [1] 7000 16970 17413
## [1] 7100 17149 17413
## [1] 7200 17328 17413
## [1] 7300 17506 17413
## [1] 7400 17684 17413
## [1] 7500 17862 17413
## [1] 7600 18119 17413
## [1] 7700 18297 17413
## [1] 7800 18476 17413
## [1] 7900 18655 17413## [1] 8000 18912 17413## [1] 8100 19090 17413
## [1] 8200 19269 17413
## [1] 8300 19526 17413
## [1] 8400 19704 17413
## [1] 8500 19881 17413
## [1] 8600 20058 17413
## [1] 8700 20312 17413
## [1] 8800 20488 17413
## [1] 8900 20664 17413
## [1] 9000 20916 17413
## [1] 9100 21091 17413
## [1] 9200 21266 17413
## [1] 9300 21441 17413
## [1] 9400 21691 17413
## [1] 9500 21866 17413
## [1] 9600 22041 17413
## [1] 9700 22218 17413
## [1] 9800 22468 17413
## [1] 9900 22643 17413
## [1] 10000 22817 17413
## [1] 10100 22992 17413
## [1] 10200 23240 17413
## [1] 10300 23414 17413
## [1] 10400 23662 17413
## [1] 10500 23836 17413
## [1] 10600 24009 17413
## [1] 10700 24256 17413
## [1] 10800 24430 17413
## [1] 10900 24679 17413
## [1] 11000 24927 17413
## [1] 11100 25250 17413
## [1] 11200 25500 17413
## [1] 11300 25823 17413
## [1] 11400 26075 17413
## [1] 11500 26397 17413
## [1] 11600 26790 17413
## [1] 11700 27106 17413
## [1] 11800 27494 17413
## [1] 11900 27959 17413
## [1] 12000 28351 17413
## [1] 12100 28892 17413
## [1] 12200 29303 17413## [1] 12300 29551 17413
## [1] 12400 29721 17413## [1] 12500 29857 17413
## [1] 12600 30024 17413
## [1] 12700 30174 17413
## [1] 12800 30333 17413
## [1] 12900 30505 17413
## [1] 13000 30643 17413
## [1] 13100 30827 17413
## [1] 13200 30971 17413
## [1] 13300 31118 17413
## [1] 13400 31311 17413
## [1] 13500 31458 17413
## [1] 13600 31658 17413
## [1] 13700 31811 17413
## [1] 13800 32021 17413
## [1] 13900 32177 17413
## [1] 14000 32392 17413
## [1] 14100 32552 17413
## [1] 14200 32775 17413
## [1] 14300 32941 17413
## [1] 14400 33112 17413
## [1] 14500 33355 17413
## [1] 14600 33534 17413
## [1] 14700 33790 17413
## [1] 14800 33968 17413
## [1] 14900 34146 17413
## [1] 15000 34400 17413
## [1] 15100 34577 17413
## [1] 15200 34753 17413
## [1] 15300 35004 17413
## [1] 15400 35179 17413
## [1] 15500 35353 17413
## [1] 15600 35602 17413
## [1] 15700 35777 17413
## [1] 15800 35952 17413
## [1] 15900 36129 17413
## [1] 16000 36377 17413
## [1] 16100 36553 17413
## [1] 16200 36728 17413
## [1] 16300 36975 17413
## [1] 16400 37148 17413
## [1] 16500 37391 17413
## [1] 16600 37635 17413
## [1] 16700 37877 17413
## [1] 16800 38188 17413
## [1] 16900 38498 17413
## [1] 17000 38808 17413
## [1] 17100 39119 17413
## [1] 17200 39428 17413
## [1] 17300 39759 17413
## [1] 17400 40107 17413NN.sigma_2100_8.5_today_S <- loop_sigma_D( B1[S_hem,], A1[S_hem,], C, "_2100_8.5_today")## [1] 100 140 22249
## [1] 200 298 22249
## [1] 300 459 22249
## [1] 400 622 22249
## [1] 500 765 22249
## [1] 600 909 22249
## [1] 700 1052 22249
## [1] 800 1195 22249
## [1] 900 1345 22249
## [1] 1000 1468 22249
## [1] 1100 1619 22249
## [1] 1200 1763 22249
## [1] 1300 1923 22249
## [1] 1400 2093 22249
## [1] 1500 2228 22249
## [1] 1600 2401 22249
## [1] 1700 2539 22249
## [1] 1800 2718 22249
## [1] 1900 2858 22249
## [1] 2000 3040 22249
## [1] 2100 3182 22249
## [1] 2200 3372 22249
## [1] 2300 3518 22249
## [1] 2400 3710 22249
## [1] 2500 3856 22249
## [1] 2600 4047 22249
## [1] 2700 4190 22249
## [1] 2800 4368 22249
## [1] 2900 4506 22249
## [1] 3000 4678 22249
## [1] 3100 4810 22249
## [1] 3200 4967 22249
## [1] 3300 5092 22249
## [1] 3400 5240 22249
## [1] 3500 5367 22249
## [1] 3600 5531 22249
## [1] 3700 5666 22249
## [1] 3800 5834 22249
## [1] 3900 5969 22249
## [1] 4000 6141 22249
## [1] 4100 6278 22249
## [1] 4200 6451 22249
## [1] 4300 6586 22249
## [1] 4400 6753 22249
## [1] 4500 6882 22249
## [1] 4600 7031 22249
## [1] 4700 7149 22249
## [1] 4800 7275 22249
## [1] 4900 7388 22249
## [1] 5000 7514 22249
## [1] 5100 7640 22249
## [1] 5200 7753 22249
## [1] 5300 7882 22249
## [1] 5400 8010 22249
## [1] 5500 8139 22249## [1] 5600 8256 22249
## [1] 5700 8404 22249
## [1] 5800 8538 22249
## [1] 5900 8680 22249
## [1] 6000 8864 22249
## [1] 6100 9073 22249
## [1] 6200 9259 22249
## [1] 6300 9501 22249
## [1] 6400 9758 22249
## [1] 6500 10019 22249
## [1] 6600 10336 22249
## [1] 6700 10614 22249
## [1] 6800 10899 22249
## [1] 6900 11128 22249
## [1] 7000 11363 22249
## [1] 7100 11532 22249
## [1] 7200 11769 22249
## [1] 7300 11938 22249
## [1] 7400 12178 22249
## [1] 7500 12350 22249
## [1] 7600 12523 22249
## [1] 7700 12772 22249
## [1] 7800 12948 22249
## [1] 7900 13203 22249
## [1] 8000 13382 22249
## [1] 8100 13561 22249
## [1] 8200 13824 22249
## [1] 8300 14012 22249
## [1] 8400 14279 22249
## [1] 8500 14462 22249
## [1] 8600 14649 22249
## [1] 8700 14915 22249
## [1] 8800 15102 22249
## [1] 8900 15289 22249
## [1] 9000 15556 22249
## [1] 9100 15740 22249
## [1] 9200 15926 22249
## [1] 9300 16112 22249
## [1] 9400 16386 22249
## [1] 9500 16573 22249
## [1] 9600 16760 22249
## [1] 9700 16947 22249
## [1] 9800 17224 22249
## [1] 9900 17414 22249
## [1] 10000 17604 22249
## [1] 10100 17881 22249
## [1] 10200 18065 22249
## [1] 10300 18247 22249
## [1] 10400 18428 22249
## [1] 10500 18688 22249
## [1] 10600 18867 22249
## [1] 10700 19044 22249
## [1] 10800 19298 22249
## [1] 10900 19475 22249
## [1] 11000 19652 22249
## [1] 11100 19828 22249
## [1] 11200 20081 22249
## [1] 11300 20258 22249
## [1] 11400 20436 22249
## [1] 11500 20694 22249
## [1] 11600 20874 22249
## [1] 11700 21054 22249
## [1] 11800 21314 22249
## [1] 11900 21494 22249
## [1] 12000 21674 22249
## [1] 12100 21934 22249
## [1] 12200 22115 22249
## [1] 12300 22295 22249
## [1] 12400 22553 22249
## [1] 12500 22731 22249
## [1] 12600 22909 22249
## [1] 12700 23161 22249
## [1] 12800 23335 22249
## [1] 12900 23508 22249
## [1] 13000 23753 22249
## [1] 13100 23923 22249
## [1] 13200 24092 22249
## [1] 13300 24326 22249
## [1] 13400 24485 22249
## [1] 13500 24690 22249
## [1] 13600 24839 22249
## [1] 13700 24987 22249
## [1] 13800 25178 22249
## [1] 13900 25322 22249
## [1] 14000 25465 22249
## [1] 14100 25646 22249
## [1] 14200 25785 22249
## [1] 14300 25923 22249
## [1] 14400 26095 22249
## [1] 14500 26228 22249
## [1] 14600 26360 22249
## [1] 14700 26522 22249
## [1] 14800 26652 22249
## [1] 14900 26811 22249
## [1] 15000 26938 22249
## [1] 15100 27065 22249
## [1] 15200 27218 22249
## [1] 15300 27344 22249
## [1] 15400 27470 22249
## [1] 15500 27619 22249
## [1] 15600 27743 22249
## [1] 15700 27890 22249
## [1] 15800 28011 22249
## [1] 15900 28132 22249
## [1] 16000 28274 22249
## [1] 16100 28394 22249
## [1] 16200 28514 22249
## [1] 16300 28648 22249
## [1] 16400 28764 22249
## [1] 16500 28898 22249
## [1] 16600 29022 22249
## [1] 16700 29172 22249
## [1] 16800 29327 22249
## [1] 16900 29498 22249
## [1] 17000 29708 22249
## [1] 17100 30060 22249
## [1] 17200 30380 22249
## [1] 17300 30680 22249
## [1] 17400 30923 22249
## [1] 17500 31234 22249
## [1] 17600 31474 22249
## [1] 17700 31716 22249
## [1] 17800 31954 22249
## [1] 17900 32192 22249
## [1] 18000 32359 22249
## [1] 18100 32598 22249
## [1] 18200 32844 22249
## [1] 18300 33019 22249
## [1] 18400 33266 22249
## [1] 18500 33439 22249
## [1] 18600 33616 22249
## [1] 18700 33863 22249
## [1] 18800 34038 22249
## [1] 18900 34213 22249
## [1] 19000 34463 22249
## [1] 19100 34638 22249
## [1] 19200 34814 22249
## [1] 19300 35066 22249
## [1] 19400 35242 22249
## [1] 19500 35419 22249
## [1] 19600 35678 22249
## [1] 19700 35859 22249
## [1] 19800 36042 22249
## [1] 19900 36312 22249
## [1] 20000 36490 22249
## [1] 20100 36664 22249
## [1] 20200 36908 22249
## [1] 20300 37078 22249
## [1] 20400 37310 22249
## [1] 20500 37464 22249
## [1] 20600 37612 22249
## [1] 20700 37801 22249
## [1] 20800 37943 22249
## [1] 20900 38122 22249
## [1] 21000 38258 22249
## [1] 21100 38428 22249
## [1] 21200 38563 22249
## [1] 21300 38697 22249
## [1] 21400 38864 22249
## [1] 21500 38997 22249
## [1] 21600 39161 22249
## [1] 21700 39292 22249
## [1] 21800 39446 22249
## [1] 21900 39572 22249
## [1] 22000 39721 22249
## [1] 22100 39866 22249
## [1] 22200 40030 22249NN.sigma_2100_8.5_today <- rbind(NN.sigma_2100_8.5_today_N, NN.sigma_2100_8.5_today_S)
final_dat7 <- full_join(NN.sigma_2100_8.5_today, final_dat6, by="No")
head(final_dat7)## No NN.sigma_2100_8.5_today NN.station_2100_8.5_today NN.Mdist_2100_8.5_today
## 1 36 8.292361 40111 10.055106
## 2 37 8.292361 40111 10.340310
## 3 38 8.292361 40111 10.400029
## 4 39 8.292361 40111 10.280321
## 5 40 8.292361 40111 10.091702
## 6 41 8.292361 40111 9.955604
## numPCs_2100_8.5_today NN.sigma_today_2100_8.5 NN.station_today_2100_8.5
## 1 2 0.3501681236 29638
## 2 2 0.1043489359 29695
## 3 2 0.0015820951 29638
## 4 2 0.0006816715 28978
## 5 2 0.1188184652 28440
## 6 2 0.5463419465 28254
## NN.Mdist_today_2100_8.5 numPCs_today_2100_8.5 lat long NN.sigma_1800_today
## 1 0.79989070 2 70.5 20.5 0.033725924
## 2 0.41656957 2 71.5 20.5 0.015960753
## 3 0.05026184 2 72.5 20.5 0.010016616
## 4 0.03298615 2 73.5 20.5 0.001600598
## 5 0.44577407 2 74.5 20.5 0.018324556
## 6 1.03579199 2 75.5 20.5 0.004808260
## NN.station_1800_today NN.Mdist_1800_today numPCs_1800_today
## 1 30751 0.23355016 2
## 2 30751 0.16010043 2
## 3 31111 0.12668130 2
## 4 14563 0.05055509 2
## 5 1958 0.17162744 2
## 6 17177 0.08767895 2
## NN.sigma_today_1800 NN.station_today_1800 NN.Mdist_today_1800
## 1 0.836709022 31931 1.3486467
## 2 0.894133247 794 1.4077485
## 3 0.661532267 861 1.1633920
## 4 0.265489647 39602 0.6854521
## 5 0.006257156 39603 0.1000496
## 6 0.036266868 39426 0.2423104
## numPCs_today_1800 lat_1800_today long_1800_today latshift_1800_today latshift
## 1 2 44.5 302.5 26 26
## 2 2 44.5 302.5 27 27
## 3 2 82.5 305.5 -10 -10
## 4 2 53.5 172.5 20 20
## 5 2 69.5 47.5 5 5
## 6 2 61.5 188.5 14 14
## lat_today_2100_8.5 long_today_2100_8.5
## 1 89.5 287.5
## 2 89.5 288.5
## 3 89.5 287.5
## 4 85.5 279.5
## 5 82.5 273.5
## 6 82.5 271.5Plot_nonInt(final_dat7$lat, final_dat7$long,
final_dat7$NN.sigma_2100_8.5_today, world, "sigma dis. 2000 in 2100 8.5")
#
# stationInfo7 <- stationInfo
# names(stationInfo7) <- paste0(names(stationInfo),"_2100_8.5_today")
# names(stationInfo7)[1] <- "NN.station_2100_8.5_today"
# head(stationInfo7)
# final_dat8 <- left_join(final_dat7, stationInfo7)
final_dat8 <- final_dat72100 RCP 4.5 against today (novelty)
#--------------------------------
### 2100 RCP 4.5 against today (novelty) ####
#--------------------------------
identical(norm_2000$No, norm_2100_4.5$No)## [1] TRUE# A1 the same as above
B2 <- norm_2100_4.5
# 2070-2100
head(B2)## No SST Arag pH
## 1 1 -0.8282522 0.04708031 7.879460
## 2 2 -0.7305017 0.04380220 7.891665
## 3 3 -0.6256332 0.04038376 7.901365
## 4 4 -0.5128039 0.03683358 7.908494
## 5 5 -0.3908187 0.03319882 7.913204
## 6 6 -0.2581154 0.02960423 7.915740dim(B2)## [1] 39662 4# same columns as A
# sanity check to make sure stations in right order
identical(A1$No, B2$No) # should be true## [1] TRUEidentical(A1$No, sort(A1$No)) # should also be true## [1] TRUE# C is the same defined above
NN.sigma_today_2100_4.5_N <- loop_sigma_D(A1[N_hem,], B2[N_hem,], C, "_today_2100_4.5")## [1] 100 317 17413
## [1] 200 610 17413
## [1] 300 934 17413
## [1] 400 1247 17413
## [1] 500 1683 17413
## [1] 600 2053 17413
## [1] 700 2356 17413
## [1] 800 2660 17413
## [1] 900 2971 17413
## [1] 1000 3213 17413
## [1] 1100 3452 17413
## [1] 1200 3759 17413
## [1] 1300 3996 17413
## [1] 1400 4236 17413
## [1] 1500 4540 17413
## [1] 1600 4847 17413
## [1] 1700 5220 17413
## [1] 1800 5530 17413
## [1] 1900 5830 17413
## [1] 2000 6129 17413
## [1] 2100 6427 17413
## [1] 2200 6726 17413
## [1] 2300 7099 17413
## [1] 2400 7694 17413
## [1] 2500 8267 17413
## [1] 2600 8629 17413
## [1] 2700 8887 17413
## [1] 2800 9069 17413
## [1] 2900 9286 17413
## [1] 3000 9454 17413
## [1] 3100 9622 17413
## [1] 3200 9790 17413
## [1] 3300 9954 17413
## [1] 3400 10119 17413
## [1] 3500 10287 17413
## [1] 3600 10458 17413
## [1] 3700 10592 17413
## [1] 3800 10764 17413
## [1] 3900 10945 17413
## [1] 4000 11095 17413
## [1] 4100 11309 17413
## [1] 4200 11469 17413
## [1] 4300 11702 17413
## [1] 4400 11871 17413
## [1] 4500 12109 17413
## [1] 4600 12278 17413
## [1] 4700 12450 17413
## [1] 4800 12692 17413
## [1] 4900 12863 17413
## [1] 5000 13035 17413
## [1] 5100 13280 17413
## [1] 5200 13453 17413
## [1] 5300 13627 17413
## [1] 5400 13878 17413
## [1] 5500 14056 17413
## [1] 5600 14235 17413
## [1] 5700 14414 17413
## [1] 5800 14670 17413
## [1] 5900 14848 17413
## [1] 6000 15026 17413
## [1] 6100 15204 17413
## [1] 6200 15384 17413
## [1] 6300 15641 17413
## [1] 6400 15819 17413
## [1] 6500 15998 17413
## [1] 6600 16176 17413
## [1] 6700 16356 17413
## [1] 6800 16534 17413
## [1] 6900 16791 17413
## [1] 7000 16970 17413
## [1] 7100 17149 17413
## [1] 7200 17328 17413
## [1] 7300 17506 17413
## [1] 7400 17684 17413
## [1] 7500 17862 17413
## [1] 7600 18119 17413
## [1] 7700 18297 17413
## [1] 7800 18476 17413
## [1] 7900 18655 17413## [1] 8000 18912 17413## [1] 8100 19090 17413
## [1] 8200 19269 17413
## [1] 8300 19526 17413
## [1] 8400 19704 17413
## [1] 8500 19881 17413
## [1] 8600 20058 17413
## [1] 8700 20312 17413
## [1] 8800 20488 17413
## [1] 8900 20664 17413
## [1] 9000 20916 17413
## [1] 9100 21091 17413
## [1] 9200 21266 17413
## [1] 9300 21441 17413
## [1] 9400 21691 17413
## [1] 9500 21866 17413
## [1] 9600 22041 17413
## [1] 9700 22218 17413
## [1] 9800 22468 17413
## [1] 9900 22643 17413
## [1] 10000 22817 17413
## [1] 10100 22992 17413
## [1] 10200 23240 17413
## [1] 10300 23414 17413
## [1] 10400 23662 17413
## [1] 10500 23836 17413
## [1] 10600 24009 17413
## [1] 10700 24256 17413
## [1] 10800 24430 17413
## [1] 10900 24679 17413
## [1] 11000 24927 17413
## [1] 11100 25250 17413
## [1] 11200 25500 17413
## [1] 11300 25823 17413
## [1] 11400 26075 17413
## [1] 11500 26397 17413
## [1] 11600 26790 17413
## [1] 11700 27106 17413
## [1] 11800 27494 17413
## [1] 11900 27959 17413
## [1] 12000 28351 17413
## [1] 12100 28892 17413
## [1] 12200 29303 17413## [1] 12300 29551 17413
## [1] 12400 29721 17413## [1] 12500 29857 17413
## [1] 12600 30024 17413
## [1] 12700 30174 17413
## [1] 12800 30333 17413
## [1] 12900 30505 17413
## [1] 13000 30643 17413
## [1] 13100 30827 17413
## [1] 13200 30971 17413
## [1] 13300 31118 17413
## [1] 13400 31311 17413
## [1] 13500 31458 17413
## [1] 13600 31658 17413
## [1] 13700 31811 17413
## [1] 13800 32021 17413
## [1] 13900 32177 17413
## [1] 14000 32392 17413
## [1] 14100 32552 17413
## [1] 14200 32775 17413
## [1] 14300 32941 17413
## [1] 14400 33112 17413
## [1] 14500 33355 17413
## [1] 14600 33534 17413
## [1] 14700 33790 17413
## [1] 14800 33968 17413
## [1] 14900 34146 17413
## [1] 15000 34400 17413
## [1] 15100 34577 17413
## [1] 15200 34753 17413
## [1] 15300 35004 17413
## [1] 15400 35179 17413
## [1] 15500 35353 17413
## [1] 15600 35602 17413
## [1] 15700 35777 17413
## [1] 15800 35952 17413
## [1] 15900 36129 17413
## [1] 16000 36377 17413
## [1] 16100 36553 17413
## [1] 16200 36728 17413
## [1] 16300 36975 17413
## [1] 16400 37148 17413
## [1] 16500 37391 17413
## [1] 16600 37635 17413
## [1] 16700 37877 17413
## [1] 16800 38188 17413
## [1] 16900 38498 17413
## [1] 17000 38808 17413
## [1] 17100 39119 17413
## [1] 17200 39428 17413
## [1] 17300 39759 17413
## [1] 17400 40107 17413NN.sigma_today_2100_4.5_S <- loop_sigma_D(A1[S_hem,], B2[S_hem,], C, "_today_2100_4.5")## [1] 100 140 22249
## [1] 200 298 22249
## [1] 300 459 22249
## [1] 400 622 22249
## [1] 500 765 22249
## [1] 600 909 22249
## [1] 700 1052 22249
## [1] 800 1195 22249
## [1] 900 1345 22249
## [1] 1000 1468 22249
## [1] 1100 1619 22249
## [1] 1200 1763 22249
## [1] 1300 1923 22249
## [1] 1400 2093 22249
## [1] 1500 2228 22249
## [1] 1600 2401 22249
## [1] 1700 2539 22249
## [1] 1800 2718 22249
## [1] 1900 2858 22249
## [1] 2000 3040 22249
## [1] 2100 3182 22249
## [1] 2200 3372 22249
## [1] 2300 3518 22249
## [1] 2400 3710 22249
## [1] 2500 3856 22249
## [1] 2600 4047 22249
## [1] 2700 4190 22249
## [1] 2800 4368 22249
## [1] 2900 4506 22249
## [1] 3000 4678 22249
## [1] 3100 4810 22249
## [1] 3200 4967 22249
## [1] 3300 5092 22249
## [1] 3400 5240 22249
## [1] 3500 5367 22249
## [1] 3600 5531 22249
## [1] 3700 5666 22249
## [1] 3800 5834 22249
## [1] 3900 5969 22249
## [1] 4000 6141 22249
## [1] 4100 6278 22249
## [1] 4200 6451 22249
## [1] 4300 6586 22249
## [1] 4400 6753 22249
## [1] 4500 6882 22249
## [1] 4600 7031 22249
## [1] 4700 7149 22249
## [1] 4800 7275 22249
## [1] 4900 7388 22249
## [1] 5000 7514 22249
## [1] 5100 7640 22249
## [1] 5200 7753 22249
## [1] 5300 7882 22249
## [1] 5400 8010 22249
## [1] 5500 8139 22249## [1] 5600 8256 22249
## [1] 5700 8404 22249
## [1] 5800 8538 22249
## [1] 5900 8680 22249
## [1] 6000 8864 22249
## [1] 6100 9073 22249
## [1] 6200 9259 22249
## [1] 6300 9501 22249
## [1] 6400 9758 22249
## [1] 6500 10019 22249
## [1] 6600 10336 22249
## [1] 6700 10614 22249
## [1] 6800 10899 22249
## [1] 6900 11128 22249
## [1] 7000 11363 22249
## [1] 7100 11532 22249
## [1] 7200 11769 22249
## [1] 7300 11938 22249
## [1] 7400 12178 22249
## [1] 7500 12350 22249
## [1] 7600 12523 22249
## [1] 7700 12772 22249
## [1] 7800 12948 22249
## [1] 7900 13203 22249
## [1] 8000 13382 22249
## [1] 8100 13561 22249
## [1] 8200 13824 22249
## [1] 8300 14012 22249
## [1] 8400 14279 22249
## [1] 8500 14462 22249
## [1] 8600 14649 22249
## [1] 8700 14915 22249
## [1] 8800 15102 22249
## [1] 8900 15289 22249
## [1] 9000 15556 22249
## [1] 9100 15740 22249
## [1] 9200 15926 22249
## [1] 9300 16112 22249
## [1] 9400 16386 22249
## [1] 9500 16573 22249
## [1] 9600 16760 22249
## [1] 9700 16947 22249
## [1] 9800 17224 22249
## [1] 9900 17414 22249
## [1] 10000 17604 22249
## [1] 10100 17881 22249
## [1] 10200 18065 22249
## [1] 10300 18247 22249
## [1] 10400 18428 22249
## [1] 10500 18688 22249
## [1] 10600 18867 22249
## [1] 10700 19044 22249
## [1] 10800 19298 22249
## [1] 10900 19475 22249
## [1] 11000 19652 22249
## [1] 11100 19828 22249
## [1] 11200 20081 22249
## [1] 11300 20258 22249
## [1] 11400 20436 22249
## [1] 11500 20694 22249
## [1] 11600 20874 22249
## [1] 11700 21054 22249
## [1] 11800 21314 22249
## [1] 11900 21494 22249
## [1] 12000 21674 22249
## [1] 12100 21934 22249
## [1] 12200 22115 22249
## [1] 12300 22295 22249
## [1] 12400 22553 22249
## [1] 12500 22731 22249
## [1] 12600 22909 22249
## [1] 12700 23161 22249
## [1] 12800 23335 22249
## [1] 12900 23508 22249
## [1] 13000 23753 22249
## [1] 13100 23923 22249
## [1] 13200 24092 22249
## [1] 13300 24326 22249
## [1] 13400 24485 22249
## [1] 13500 24690 22249
## [1] 13600 24839 22249
## [1] 13700 24987 22249
## [1] 13800 25178 22249
## [1] 13900 25322 22249
## [1] 14000 25465 22249
## [1] 14100 25646 22249
## [1] 14200 25785 22249
## [1] 14300 25923 22249
## [1] 14400 26095 22249
## [1] 14500 26228 22249
## [1] 14600 26360 22249
## [1] 14700 26522 22249
## [1] 14800 26652 22249
## [1] 14900 26811 22249
## [1] 15000 26938 22249
## [1] 15100 27065 22249
## [1] 15200 27218 22249
## [1] 15300 27344 22249
## [1] 15400 27470 22249
## [1] 15500 27619 22249
## [1] 15600 27743 22249
## [1] 15700 27890 22249
## [1] 15800 28011 22249
## [1] 15900 28132 22249
## [1] 16000 28274 22249
## [1] 16100 28394 22249
## [1] 16200 28514 22249
## [1] 16300 28648 22249
## [1] 16400 28764 22249
## [1] 16500 28898 22249
## [1] 16600 29022 22249
## [1] 16700 29172 22249
## [1] 16800 29327 22249
## [1] 16900 29498 22249
## [1] 17000 29708 22249
## [1] 17100 30060 22249
## [1] 17200 30380 22249
## [1] 17300 30680 22249
## [1] 17400 30923 22249
## [1] 17500 31234 22249
## [1] 17600 31474 22249
## [1] 17700 31716 22249
## [1] 17800 31954 22249
## [1] 17900 32192 22249
## [1] 18000 32359 22249
## [1] 18100 32598 22249
## [1] 18200 32844 22249
## [1] 18300 33019 22249
## [1] 18400 33266 22249
## [1] 18500 33439 22249
## [1] 18600 33616 22249
## [1] 18700 33863 22249
## [1] 18800 34038 22249
## [1] 18900 34213 22249
## [1] 19000 34463 22249
## [1] 19100 34638 22249
## [1] 19200 34814 22249
## [1] 19300 35066 22249
## [1] 19400 35242 22249
## [1] 19500 35419 22249
## [1] 19600 35678 22249
## [1] 19700 35859 22249
## [1] 19800 36042 22249
## [1] 19900 36312 22249
## [1] 20000 36490 22249
## [1] 20100 36664 22249
## [1] 20200 36908 22249
## [1] 20300 37078 22249
## [1] 20400 37310 22249
## [1] 20500 37464 22249
## [1] 20600 37612 22249
## [1] 20700 37801 22249
## [1] 20800 37943 22249
## [1] 20900 38122 22249
## [1] 21000 38258 22249
## [1] 21100 38428 22249
## [1] 21200 38563 22249
## [1] 21300 38697 22249
## [1] 21400 38864 22249
## [1] 21500 38997 22249
## [1] 21600 39161 22249
## [1] 21700 39292 22249
## [1] 21800 39446 22249
## [1] 21900 39572 22249
## [1] 22000 39721 22249
## [1] 22100 39866 22249
## [1] 22200 40030 22249NN.sigma_today_2100_4.5 <- rbind(NN.sigma_today_2100_4.5_N, NN.sigma_today_2100_4.5_S)
which(is.infinite(NN.sigma_today_2100_4.5$NN.sigma_today_2100_4.5))## integer(0)which(is.na(NN.sigma_today_2100_4.5$NN.sigma_today_2100_4.5))## integer(0)final_dat9 <- full_join(NN.sigma_today_2100_4.5, final_dat8, by="No")
head(final_dat9)## No NN.sigma_today_2100_4.5 NN.station_today_2100_4.5 NN.Mdist_today_2100_4.5
## 1 36 0.260091777 31342 0.67775158
## 2 37 0.150575943 31338 0.50493740
## 3 38 0.005534165 31456 0.09407847
## 4 39 0.139495959 30347 0.48495881
## 5 40 0.011775140 3437 0.13740012
## 6 41 0.098007656 4227 0.40321232
## numPCs_today_2100_4.5 NN.sigma_2100_8.5_today NN.station_2100_8.5_today
## 1 2 8.292361 40111
## 2 2 8.292361 40111
## 3 2 8.292361 40111
## 4 2 8.292361 40111
## 5 2 8.292361 40111
## 6 2 8.292361 40111
## NN.Mdist_2100_8.5_today numPCs_2100_8.5_today NN.sigma_today_2100_8.5
## 1 10.055106 2 0.3501681236
## 2 10.340310 2 0.1043489359
## 3 10.400029 2 0.0015820951
## 4 10.280321 2 0.0006816715
## 5 10.091702 2 0.1188184652
## 6 9.955604 2 0.5463419465
## NN.station_today_2100_8.5 NN.Mdist_today_2100_8.5 numPCs_today_2100_8.5 lat
## 1 29638 0.79989070 2 70.5
## 2 29695 0.41656957 2 71.5
## 3 29638 0.05026184 2 72.5
## 4 28978 0.03298615 2 73.5
## 5 28440 0.44577407 2 74.5
## 6 28254 1.03579199 2 75.5
## long NN.sigma_1800_today NN.station_1800_today NN.Mdist_1800_today
## 1 20.5 0.033725924 30751 0.23355016
## 2 20.5 0.015960753 30751 0.16010043
## 3 20.5 0.010016616 31111 0.12668130
## 4 20.5 0.001600598 14563 0.05055509
## 5 20.5 0.018324556 1958 0.17162744
## 6 20.5 0.004808260 17177 0.08767895
## numPCs_1800_today NN.sigma_today_1800 NN.station_today_1800
## 1 2 0.836709022 31931
## 2 2 0.894133247 794
## 3 2 0.661532267 861
## 4 2 0.265489647 39602
## 5 2 0.006257156 39603
## 6 2 0.036266868 39426
## NN.Mdist_today_1800 numPCs_today_1800 lat_1800_today long_1800_today
## 1 1.3486467 2 44.5 302.5
## 2 1.4077485 2 44.5 302.5
## 3 1.1633920 2 82.5 305.5
## 4 0.6854521 2 53.5 172.5
## 5 0.1000496 2 69.5 47.5
## 6 0.2423104 2 61.5 188.5
## latshift_1800_today latshift lat_today_2100_8.5 long_today_2100_8.5
## 1 26 26 89.5 287.5
## 2 27 27 89.5 288.5
## 3 -10 -10 89.5 287.5
## 4 20 20 85.5 279.5
## 5 5 5 82.5 273.5
## 6 14 14 82.5 271.5# Visualize
Plot_nonInt(final_dat9$lat, final_dat9$long,
final_dat9$NN.sigma_today_2100_4.5, world, "sig nov 2100 4.5")
stationInfo9 <- stationInfo
names(stationInfo9) <- paste0(names(stationInfo),"_today_2100_4.5")
names(stationInfo9)[1] <- "NN.station_today_2100_4.5"
head(stationInfo9)## NN.station_today_2100_4.5 lat_today_2100_4.5 long_today_2100_4.5
## 1 1 -69.5 20.5
## 2 2 -68.5 20.5
## 3 3 -67.5 20.5
## 4 4 -66.5 20.5
## 5 5 -65.5 20.5
## 6 6 -64.5 20.5final_dat10 <- left_join(final_dat9, stationInfo9)## Joining, by = "NN.station_today_2100_4.5"dim(final_dat10)## [1] 39662 31head(final_dat10)## No NN.sigma_today_2100_4.5 NN.station_today_2100_4.5 NN.Mdist_today_2100_4.5
## 1 36 0.260091777 31342 0.67775158
## 2 37 0.150575943 31338 0.50493740
## 3 38 0.005534165 31456 0.09407847
## 4 39 0.139495959 30347 0.48495881
## 5 40 0.011775140 3437 0.13740012
## 6 41 0.098007656 4227 0.40321232
## numPCs_today_2100_4.5 NN.sigma_2100_8.5_today NN.station_2100_8.5_today
## 1 2 8.292361 40111
## 2 2 8.292361 40111
## 3 2 8.292361 40111
## 4 2 8.292361 40111
## 5 2 8.292361 40111
## 6 2 8.292361 40111
## NN.Mdist_2100_8.5_today numPCs_2100_8.5_today NN.sigma_today_2100_8.5
## 1 10.055106 2 0.3501681236
## 2 10.340310 2 0.1043489359
## 3 10.400029 2 0.0015820951
## 4 10.280321 2 0.0006816715
## 5 10.091702 2 0.1188184652
## 6 9.955604 2 0.5463419465
## NN.station_today_2100_8.5 NN.Mdist_today_2100_8.5 numPCs_today_2100_8.5 lat
## 1 29638 0.79989070 2 70.5
## 2 29695 0.41656957 2 71.5
## 3 29638 0.05026184 2 72.5
## 4 28978 0.03298615 2 73.5
## 5 28440 0.44577407 2 74.5
## 6 28254 1.03579199 2 75.5
## long NN.sigma_1800_today NN.station_1800_today NN.Mdist_1800_today
## 1 20.5 0.033725924 30751 0.23355016
## 2 20.5 0.015960753 30751 0.16010043
## 3 20.5 0.010016616 31111 0.12668130
## 4 20.5 0.001600598 14563 0.05055509
## 5 20.5 0.018324556 1958 0.17162744
## 6 20.5 0.004808260 17177 0.08767895
## numPCs_1800_today NN.sigma_today_1800 NN.station_today_1800
## 1 2 0.836709022 31931
## 2 2 0.894133247 794
## 3 2 0.661532267 861
## 4 2 0.265489647 39602
## 5 2 0.006257156 39603
## 6 2 0.036266868 39426
## NN.Mdist_today_1800 numPCs_today_1800 lat_1800_today long_1800_today
## 1 1.3486467 2 44.5 302.5
## 2 1.4077485 2 44.5 302.5
## 3 1.1633920 2 82.5 305.5
## 4 0.6854521 2 53.5 172.5
## 5 0.1000496 2 69.5 47.5
## 6 0.2423104 2 61.5 188.5
## latshift_1800_today latshift lat_today_2100_8.5 long_today_2100_8.5
## 1 26 26 89.5 287.5
## 2 27 27 89.5 288.5
## 3 -10 -10 89.5 287.5
## 4 20 20 85.5 279.5
## 5 5 5 82.5 273.5
## 6 14 14 82.5 271.5
## lat_today_2100_4.5 long_today_2100_4.5
## 1 86.5 307.5
## 2 82.5 307.5
## 3 85.5 308.5
## 4 69.5 298.5
## 5 70.5 61.5
## 6 72.5 68.5today against 2100 RCP 4.5 (disappearing climate)
#--------------------------------
### today against 2100 RCP 4.5 (disappearing climate) ####
#--------------------------------
NN.sigma_2100_4.5_today_N <- loop_sigma_D(B2[N_hem,], A1[N_hem,], C, "_2100_4.5_today")## [1] 100 317 17413
## [1] 200 610 17413
## [1] 300 934 17413
## [1] 400 1247 17413
## [1] 500 1683 17413
## [1] 600 2053 17413
## [1] 700 2356 17413
## [1] 800 2660 17413
## [1] 900 2971 17413
## [1] 1000 3213 17413
## [1] 1100 3452 17413
## [1] 1200 3759 17413
## [1] 1300 3996 17413
## [1] 1400 4236 17413
## [1] 1500 4540 17413
## [1] 1600 4847 17413
## [1] 1700 5220 17413
## [1] 1800 5530 17413
## [1] 1900 5830 17413
## [1] 2000 6129 17413
## [1] 2100 6427 17413
## [1] 2200 6726 17413
## [1] 2300 7099 17413
## [1] 2400 7694 17413
## [1] 2500 8267 17413
## [1] 2600 8629 17413
## [1] 2700 8887 17413
## [1] 2800 9069 17413
## [1] 2900 9286 17413
## [1] 3000 9454 17413
## [1] 3100 9622 17413
## [1] 3200 9790 17413
## [1] 3300 9954 17413
## [1] 3400 10119 17413
## [1] 3500 10287 17413
## [1] 3600 10458 17413
## [1] 3700 10592 17413
## [1] 3800 10764 17413
## [1] 3900 10945 17413
## [1] 4000 11095 17413
## [1] 4100 11309 17413
## [1] 4200 11469 17413
## [1] 4300 11702 17413
## [1] 4400 11871 17413
## [1] 4500 12109 17413
## [1] 4600 12278 17413
## [1] 4700 12450 17413
## [1] 4800 12692 17413
## [1] 4900 12863 17413
## [1] 5000 13035 17413
## [1] 5100 13280 17413
## [1] 5200 13453 17413
## [1] 5300 13627 17413
## [1] 5400 13878 17413
## [1] 5500 14056 17413
## [1] 5600 14235 17413
## [1] 5700 14414 17413
## [1] 5800 14670 17413
## [1] 5900 14848 17413
## [1] 6000 15026 17413
## [1] 6100 15204 17413
## [1] 6200 15384 17413
## [1] 6300 15641 17413
## [1] 6400 15819 17413
## [1] 6500 15998 17413
## [1] 6600 16176 17413
## [1] 6700 16356 17413
## [1] 6800 16534 17413
## [1] 6900 16791 17413
## [1] 7000 16970 17413
## [1] 7100 17149 17413
## [1] 7200 17328 17413
## [1] 7300 17506 17413
## [1] 7400 17684 17413
## [1] 7500 17862 17413
## [1] 7600 18119 17413
## [1] 7700 18297 17413
## [1] 7800 18476 17413
## [1] 7900 18655 17413## [1] 8000 18912 17413## [1] 8100 19090 17413
## [1] 8200 19269 17413
## [1] 8300 19526 17413
## [1] 8400 19704 17413
## [1] 8500 19881 17413
## [1] 8600 20058 17413
## [1] 8700 20312 17413
## [1] 8800 20488 17413
## [1] 8900 20664 17413
## [1] 9000 20916 17413
## [1] 9100 21091 17413
## [1] 9200 21266 17413
## [1] 9300 21441 17413
## [1] 9400 21691 17413
## [1] 9500 21866 17413
## [1] 9600 22041 17413
## [1] 9700 22218 17413
## [1] 9800 22468 17413
## [1] 9900 22643 17413
## [1] 10000 22817 17413
## [1] 10100 22992 17413
## [1] 10200 23240 17413
## [1] 10300 23414 17413
## [1] 10400 23662 17413
## [1] 10500 23836 17413
## [1] 10600 24009 17413
## [1] 10700 24256 17413
## [1] 10800 24430 17413
## [1] 10900 24679 17413
## [1] 11000 24927 17413
## [1] 11100 25250 17413
## [1] 11200 25500 17413
## [1] 11300 25823 17413
## [1] 11400 26075 17413
## [1] 11500 26397 17413
## [1] 11600 26790 17413
## [1] 11700 27106 17413
## [1] 11800 27494 17413
## [1] 11900 27959 17413
## [1] 12000 28351 17413
## [1] 12100 28892 17413
## [1] 12200 29303 17413## [1] 12300 29551 17413
## [1] 12400 29721 17413## [1] 12500 29857 17413
## [1] 12600 30024 17413
## [1] 12700 30174 17413
## [1] 12800 30333 17413
## [1] 12900 30505 17413
## [1] 13000 30643 17413
## [1] 13100 30827 17413
## [1] 13200 30971 17413
## [1] 13300 31118 17413
## [1] 13400 31311 17413
## [1] 13500 31458 17413
## [1] 13600 31658 17413
## [1] 13700 31811 17413
## [1] 13800 32021 17413
## [1] 13900 32177 17413
## [1] 14000 32392 17413
## [1] 14100 32552 17413
## [1] 14200 32775 17413
## [1] 14300 32941 17413
## [1] 14400 33112 17413
## [1] 14500 33355 17413
## [1] 14600 33534 17413
## [1] 14700 33790 17413
## [1] 14800 33968 17413
## [1] 14900 34146 17413
## [1] 15000 34400 17413
## [1] 15100 34577 17413
## [1] 15200 34753 17413
## [1] 15300 35004 17413
## [1] 15400 35179 17413
## [1] 15500 35353 17413
## [1] 15600 35602 17413
## [1] 15700 35777 17413
## [1] 15800 35952 17413
## [1] 15900 36129 17413
## [1] 16000 36377 17413
## [1] 16100 36553 17413
## [1] 16200 36728 17413
## [1] 16300 36975 17413
## [1] 16400 37148 17413
## [1] 16500 37391 17413
## [1] 16600 37635 17413
## [1] 16700 37877 17413
## [1] 16800 38188 17413
## [1] 16900 38498 17413
## [1] 17000 38808 17413
## [1] 17100 39119 17413
## [1] 17200 39428 17413
## [1] 17300 39759 17413
## [1] 17400 40107 17413NN.sigma_2100_4.5_today_S <- loop_sigma_D(B2[S_hem,], A1[S_hem,], C, "_2100_4.5_today")## [1] 100 140 22249
## [1] 200 298 22249
## [1] 300 459 22249
## [1] 400 622 22249
## [1] 500 765 22249
## [1] 600 909 22249
## [1] 700 1052 22249
## [1] 800 1195 22249
## [1] 900 1345 22249
## [1] 1000 1468 22249
## [1] 1100 1619 22249
## [1] 1200 1763 22249
## [1] 1300 1923 22249
## [1] 1400 2093 22249
## [1] 1500 2228 22249
## [1] 1600 2401 22249
## [1] 1700 2539 22249
## [1] 1800 2718 22249
## [1] 1900 2858 22249
## [1] 2000 3040 22249
## [1] 2100 3182 22249
## [1] 2200 3372 22249
## [1] 2300 3518 22249
## [1] 2400 3710 22249
## [1] 2500 3856 22249
## [1] 2600 4047 22249
## [1] 2700 4190 22249
## [1] 2800 4368 22249
## [1] 2900 4506 22249
## [1] 3000 4678 22249
## [1] 3100 4810 22249
## [1] 3200 4967 22249
## [1] 3300 5092 22249
## [1] 3400 5240 22249
## [1] 3500 5367 22249
## [1] 3600 5531 22249
## [1] 3700 5666 22249
## [1] 3800 5834 22249
## [1] 3900 5969 22249
## [1] 4000 6141 22249
## [1] 4100 6278 22249
## [1] 4200 6451 22249
## [1] 4300 6586 22249
## [1] 4400 6753 22249
## [1] 4500 6882 22249
## [1] 4600 7031 22249
## [1] 4700 7149 22249
## [1] 4800 7275 22249
## [1] 4900 7388 22249
## [1] 5000 7514 22249
## [1] 5100 7640 22249
## [1] 5200 7753 22249
## [1] 5300 7882 22249
## [1] 5400 8010 22249
## [1] 5500 8139 22249## [1] 5600 8256 22249
## [1] 5700 8404 22249
## [1] 5800 8538 22249
## [1] 5900 8680 22249
## [1] 6000 8864 22249
## [1] 6100 9073 22249
## [1] 6200 9259 22249
## [1] 6300 9501 22249
## [1] 6400 9758 22249
## [1] 6500 10019 22249
## [1] 6600 10336 22249
## [1] 6700 10614 22249
## [1] 6800 10899 22249
## [1] 6900 11128 22249
## [1] 7000 11363 22249
## [1] 7100 11532 22249
## [1] 7200 11769 22249
## [1] 7300 11938 22249
## [1] 7400 12178 22249
## [1] 7500 12350 22249
## [1] 7600 12523 22249
## [1] 7700 12772 22249
## [1] 7800 12948 22249
## [1] 7900 13203 22249
## [1] 8000 13382 22249
## [1] 8100 13561 22249
## [1] 8200 13824 22249
## [1] 8300 14012 22249
## [1] 8400 14279 22249
## [1] 8500 14462 22249
## [1] 8600 14649 22249
## [1] 8700 14915 22249
## [1] 8800 15102 22249
## [1] 8900 15289 22249
## [1] 9000 15556 22249
## [1] 9100 15740 22249
## [1] 9200 15926 22249
## [1] 9300 16112 22249
## [1] 9400 16386 22249
## [1] 9500 16573 22249
## [1] 9600 16760 22249
## [1] 9700 16947 22249
## [1] 9800 17224 22249
## [1] 9900 17414 22249
## [1] 10000 17604 22249
## [1] 10100 17881 22249
## [1] 10200 18065 22249
## [1] 10300 18247 22249
## [1] 10400 18428 22249
## [1] 10500 18688 22249
## [1] 10600 18867 22249
## [1] 10700 19044 22249
## [1] 10800 19298 22249
## [1] 10900 19475 22249
## [1] 11000 19652 22249
## [1] 11100 19828 22249
## [1] 11200 20081 22249
## [1] 11300 20258 22249
## [1] 11400 20436 22249
## [1] 11500 20694 22249
## [1] 11600 20874 22249
## [1] 11700 21054 22249
## [1] 11800 21314 22249
## [1] 11900 21494 22249
## [1] 12000 21674 22249
## [1] 12100 21934 22249
## [1] 12200 22115 22249
## [1] 12300 22295 22249
## [1] 12400 22553 22249
## [1] 12500 22731 22249
## [1] 12600 22909 22249
## [1] 12700 23161 22249
## [1] 12800 23335 22249
## [1] 12900 23508 22249
## [1] 13000 23753 22249
## [1] 13100 23923 22249
## [1] 13200 24092 22249
## [1] 13300 24326 22249
## [1] 13400 24485 22249
## [1] 13500 24690 22249
## [1] 13600 24839 22249
## [1] 13700 24987 22249
## [1] 13800 25178 22249
## [1] 13900 25322 22249
## [1] 14000 25465 22249
## [1] 14100 25646 22249
## [1] 14200 25785 22249
## [1] 14300 25923 22249
## [1] 14400 26095 22249
## [1] 14500 26228 22249
## [1] 14600 26360 22249
## [1] 14700 26522 22249
## [1] 14800 26652 22249
## [1] 14900 26811 22249
## [1] 15000 26938 22249
## [1] 15100 27065 22249
## [1] 15200 27218 22249
## [1] 15300 27344 22249
## [1] 15400 27470 22249
## [1] 15500 27619 22249
## [1] 15600 27743 22249
## [1] 15700 27890 22249
## [1] 15800 28011 22249
## [1] 15900 28132 22249
## [1] 16000 28274 22249
## [1] 16100 28394 22249
## [1] 16200 28514 22249
## [1] 16300 28648 22249
## [1] 16400 28764 22249
## [1] 16500 28898 22249
## [1] 16600 29022 22249
## [1] 16700 29172 22249
## [1] 16800 29327 22249
## [1] 16900 29498 22249
## [1] 17000 29708 22249
## [1] 17100 30060 22249
## [1] 17200 30380 22249
## [1] 17300 30680 22249
## [1] 17400 30923 22249
## [1] 17500 31234 22249
## [1] 17600 31474 22249
## [1] 17700 31716 22249
## [1] 17800 31954 22249
## [1] 17900 32192 22249
## [1] 18000 32359 22249
## [1] 18100 32598 22249
## [1] 18200 32844 22249
## [1] 18300 33019 22249
## [1] 18400 33266 22249
## [1] 18500 33439 22249
## [1] 18600 33616 22249
## [1] 18700 33863 22249
## [1] 18800 34038 22249
## [1] 18900 34213 22249
## [1] 19000 34463 22249
## [1] 19100 34638 22249
## [1] 19200 34814 22249
## [1] 19300 35066 22249
## [1] 19400 35242 22249
## [1] 19500 35419 22249
## [1] 19600 35678 22249
## [1] 19700 35859 22249
## [1] 19800 36042 22249
## [1] 19900 36312 22249
## [1] 20000 36490 22249
## [1] 20100 36664 22249
## [1] 20200 36908 22249
## [1] 20300 37078 22249
## [1] 20400 37310 22249
## [1] 20500 37464 22249
## [1] 20600 37612 22249
## [1] 20700 37801 22249
## [1] 20800 37943 22249
## [1] 20900 38122 22249
## [1] 21000 38258 22249
## [1] 21100 38428 22249
## [1] 21200 38563 22249
## [1] 21300 38697 22249
## [1] 21400 38864 22249
## [1] 21500 38997 22249
## [1] 21600 39161 22249
## [1] 21700 39292 22249
## [1] 21800 39446 22249
## [1] 21900 39572 22249
## [1] 22000 39721 22249
## [1] 22100 39866 22249
## [1] 22200 40030 22249NN.sigma_2100_4.5_today <- rbind(NN.sigma_2100_4.5_today_N, NN.sigma_2100_4.5_today_S)
which(is.infinite(NN.sigma_2100_4.5_today$NN.sigma_2100_4.5_today))## integer(0)which(is.na(NN.sigma_2100_4.5_today$NN.sigma_2100_4.5_today))## integer(0)final_dat11 <- full_join(NN.sigma_2100_4.5_today, final_dat10, by="No")
head(final_dat11)## No NN.sigma_2100_4.5_today NN.station_2100_4.5_today NN.Mdist_2100_4.5_today
## 1 36 3.551980 40055 3.967162
## 2 37 3.505984 40055 3.923107
## 3 38 3.315848 40055 3.741139
## 4 39 2.951652 40055 3.393209
## 5 40 2.553620 40111 3.013691
## 6 41 2.223800 40111 2.699426
## numPCs_2100_4.5_today NN.sigma_today_2100_4.5 NN.station_today_2100_4.5
## 1 2 0.260091777 31342
## 2 2 0.150575943 31338
## 3 2 0.005534165 31456
## 4 2 0.139495959 30347
## 5 2 0.011775140 3437
## 6 2 0.098007656 4227
## NN.Mdist_today_2100_4.5 numPCs_today_2100_4.5 NN.sigma_2100_8.5_today
## 1 0.67775158 2 8.292361
## 2 0.50493740 2 8.292361
## 3 0.09407847 2 8.292361
## 4 0.48495881 2 8.292361
## 5 0.13740012 2 8.292361
## 6 0.40321232 2 8.292361
## NN.station_2100_8.5_today NN.Mdist_2100_8.5_today numPCs_2100_8.5_today
## 1 40111 10.055106 2
## 2 40111 10.340310 2
## 3 40111 10.400029 2
## 4 40111 10.280321 2
## 5 40111 10.091702 2
## 6 40111 9.955604 2
## NN.sigma_today_2100_8.5 NN.station_today_2100_8.5 NN.Mdist_today_2100_8.5
## 1 0.3501681236 29638 0.79989070
## 2 0.1043489359 29695 0.41656957
## 3 0.0015820951 29638 0.05026184
## 4 0.0006816715 28978 0.03298615
## 5 0.1188184652 28440 0.44577407
## 6 0.5463419465 28254 1.03579199
## numPCs_today_2100_8.5 lat long NN.sigma_1800_today NN.station_1800_today
## 1 2 70.5 20.5 0.033725924 30751
## 2 2 71.5 20.5 0.015960753 30751
## 3 2 72.5 20.5 0.010016616 31111
## 4 2 73.5 20.5 0.001600598 14563
## 5 2 74.5 20.5 0.018324556 1958
## 6 2 75.5 20.5 0.004808260 17177
## NN.Mdist_1800_today numPCs_1800_today NN.sigma_today_1800
## 1 0.23355016 2 0.836709022
## 2 0.16010043 2 0.894133247
## 3 0.12668130 2 0.661532267
## 4 0.05055509 2 0.265489647
## 5 0.17162744 2 0.006257156
## 6 0.08767895 2 0.036266868
## NN.station_today_1800 NN.Mdist_today_1800 numPCs_today_1800 lat_1800_today
## 1 31931 1.3486467 2 44.5
## 2 794 1.4077485 2 44.5
## 3 861 1.1633920 2 82.5
## 4 39602 0.6854521 2 53.5
## 5 39603 0.1000496 2 69.5
## 6 39426 0.2423104 2 61.5
## long_1800_today latshift_1800_today latshift lat_today_2100_8.5
## 1 302.5 26 26 89.5
## 2 302.5 27 27 89.5
## 3 305.5 -10 -10 89.5
## 4 172.5 20 20 85.5
## 5 47.5 5 5 82.5
## 6 188.5 14 14 82.5
## long_today_2100_8.5 lat_today_2100_4.5 long_today_2100_4.5
## 1 287.5 86.5 307.5
## 2 288.5 82.5 307.5
## 3 287.5 85.5 308.5
## 4 279.5 69.5 298.5
## 5 273.5 70.5 61.5
## 6 271.5 72.5 68.5# Visualize
Plot_nonInt(final_dat11$lat, final_dat11$long,
final_dat11$NN.sigma_2100_4.5_today, world, "sigma dis. 2000 in 2100 4.5")
# stationInfo11 <- stationInfo
# names(stationInfo11) <- paste0(names(stationInfo),"_2100_4.5_today")
# names(stationInfo11)[1] <- "NN.station_2100_4.5_today"
# head(stationInfo11)
# final_dat12 <- left_join(final_dat11, stationInfo11)
final_dat12 <- final_dat11
dim(final_dat12)## [1] 39662 35head(final_dat12)## No NN.sigma_2100_4.5_today NN.station_2100_4.5_today NN.Mdist_2100_4.5_today
## 1 36 3.551980 40055 3.967162
## 2 37 3.505984 40055 3.923107
## 3 38 3.315848 40055 3.741139
## 4 39 2.951652 40055 3.393209
## 5 40 2.553620 40111 3.013691
## 6 41 2.223800 40111 2.699426
## numPCs_2100_4.5_today NN.sigma_today_2100_4.5 NN.station_today_2100_4.5
## 1 2 0.260091777 31342
## 2 2 0.150575943 31338
## 3 2 0.005534165 31456
## 4 2 0.139495959 30347
## 5 2 0.011775140 3437
## 6 2 0.098007656 4227
## NN.Mdist_today_2100_4.5 numPCs_today_2100_4.5 NN.sigma_2100_8.5_today
## 1 0.67775158 2 8.292361
## 2 0.50493740 2 8.292361
## 3 0.09407847 2 8.292361
## 4 0.48495881 2 8.292361
## 5 0.13740012 2 8.292361
## 6 0.40321232 2 8.292361
## NN.station_2100_8.5_today NN.Mdist_2100_8.5_today numPCs_2100_8.5_today
## 1 40111 10.055106 2
## 2 40111 10.340310 2
## 3 40111 10.400029 2
## 4 40111 10.280321 2
## 5 40111 10.091702 2
## 6 40111 9.955604 2
## NN.sigma_today_2100_8.5 NN.station_today_2100_8.5 NN.Mdist_today_2100_8.5
## 1 0.3501681236 29638 0.79989070
## 2 0.1043489359 29695 0.41656957
## 3 0.0015820951 29638 0.05026184
## 4 0.0006816715 28978 0.03298615
## 5 0.1188184652 28440 0.44577407
## 6 0.5463419465 28254 1.03579199
## numPCs_today_2100_8.5 lat long NN.sigma_1800_today NN.station_1800_today
## 1 2 70.5 20.5 0.033725924 30751
## 2 2 71.5 20.5 0.015960753 30751
## 3 2 72.5 20.5 0.010016616 31111
## 4 2 73.5 20.5 0.001600598 14563
## 5 2 74.5 20.5 0.018324556 1958
## 6 2 75.5 20.5 0.004808260 17177
## NN.Mdist_1800_today numPCs_1800_today NN.sigma_today_1800
## 1 0.23355016 2 0.836709022
## 2 0.16010043 2 0.894133247
## 3 0.12668130 2 0.661532267
## 4 0.05055509 2 0.265489647
## 5 0.17162744 2 0.006257156
## 6 0.08767895 2 0.036266868
## NN.station_today_1800 NN.Mdist_today_1800 numPCs_today_1800 lat_1800_today
## 1 31931 1.3486467 2 44.5
## 2 794 1.4077485 2 44.5
## 3 861 1.1633920 2 82.5
## 4 39602 0.6854521 2 53.5
## 5 39603 0.1000496 2 69.5
## 6 39426 0.2423104 2 61.5
## long_1800_today latshift_1800_today latshift lat_today_2100_8.5
## 1 302.5 26 26 89.5
## 2 302.5 27 27 89.5
## 3 305.5 -10 -10 89.5
## 4 172.5 20 20 85.5
## 5 47.5 5 5 82.5
## 6 188.5 14 14 82.5
## long_today_2100_8.5 lat_today_2100_4.5 long_today_2100_4.5
## 1 287.5 86.5 307.5
## 2 288.5 82.5 307.5
## 3 287.5 85.5 308.5
## 4 279.5 69.5 298.5
## 5 273.5 70.5 61.5
## 6 271.5 72.5 68.5sum(!complete.cases(final_dat12))## [1] 22249### Final dataframe ####
dat_Nov <- final_dat12
### Write to file ####
write.csv(final_dat12, paste0(results_dir, "SigmaD_Hem.csv"), row.names = FALSE)
head(dat_Nov)## No NN.sigma_2100_4.5_today NN.station_2100_4.5_today NN.Mdist_2100_4.5_today
## 1 36 3.551980 40055 3.967162
## 2 37 3.505984 40055 3.923107
## 3 38 3.315848 40055 3.741139
## 4 39 2.951652 40055 3.393209
## 5 40 2.553620 40111 3.013691
## 6 41 2.223800 40111 2.699426
## numPCs_2100_4.5_today NN.sigma_today_2100_4.5 NN.station_today_2100_4.5
## 1 2 0.260091777 31342
## 2 2 0.150575943 31338
## 3 2 0.005534165 31456
## 4 2 0.139495959 30347
## 5 2 0.011775140 3437
## 6 2 0.098007656 4227
## NN.Mdist_today_2100_4.5 numPCs_today_2100_4.5 NN.sigma_2100_8.5_today
## 1 0.67775158 2 8.292361
## 2 0.50493740 2 8.292361
## 3 0.09407847 2 8.292361
## 4 0.48495881 2 8.292361
## 5 0.13740012 2 8.292361
## 6 0.40321232 2 8.292361
## NN.station_2100_8.5_today NN.Mdist_2100_8.5_today numPCs_2100_8.5_today
## 1 40111 10.055106 2
## 2 40111 10.340310 2
## 3 40111 10.400029 2
## 4 40111 10.280321 2
## 5 40111 10.091702 2
## 6 40111 9.955604 2
## NN.sigma_today_2100_8.5 NN.station_today_2100_8.5 NN.Mdist_today_2100_8.5
## 1 0.3501681236 29638 0.79989070
## 2 0.1043489359 29695 0.41656957
## 3 0.0015820951 29638 0.05026184
## 4 0.0006816715 28978 0.03298615
## 5 0.1188184652 28440 0.44577407
## 6 0.5463419465 28254 1.03579199
## numPCs_today_2100_8.5 lat long NN.sigma_1800_today NN.station_1800_today
## 1 2 70.5 20.5 0.033725924 30751
## 2 2 71.5 20.5 0.015960753 30751
## 3 2 72.5 20.5 0.010016616 31111
## 4 2 73.5 20.5 0.001600598 14563
## 5 2 74.5 20.5 0.018324556 1958
## 6 2 75.5 20.5 0.004808260 17177
## NN.Mdist_1800_today numPCs_1800_today NN.sigma_today_1800
## 1 0.23355016 2 0.836709022
## 2 0.16010043 2 0.894133247
## 3 0.12668130 2 0.661532267
## 4 0.05055509 2 0.265489647
## 5 0.17162744 2 0.006257156
## 6 0.08767895 2 0.036266868
## NN.station_today_1800 NN.Mdist_today_1800 numPCs_today_1800 lat_1800_today
## 1 31931 1.3486467 2 44.5
## 2 794 1.4077485 2 44.5
## 3 861 1.1633920 2 82.5
## 4 39602 0.6854521 2 53.5
## 5 39603 0.1000496 2 69.5
## 6 39426 0.2423104 2 61.5
## long_1800_today latshift_1800_today latshift lat_today_2100_8.5
## 1 302.5 26 26 89.5
## 2 302.5 27 27 89.5
## 3 305.5 -10 -10 89.5
## 4 172.5 20 20 85.5
## 5 47.5 5 5 82.5
## 6 188.5 14 14 82.5
## long_today_2100_8.5 lat_today_2100_4.5 long_today_2100_4.5
## 1 287.5 86.5 307.5
## 2 288.5 82.5 307.5
## 3 287.5 85.5 308.5
## 4 279.5 69.5 298.5
## 5 273.5 70.5 61.5
## 6 271.5 72.5 68.5for_liqing <- data.frame(No=dat_Nov$No, lat=dat_Nov$lat,long=dat_Nov$long,
A=dat_Nov$NN.sigma_today_1800,
B=dat_Nov$NN.sigma_1800_today,
C=dat_Nov$NN.sigma_2100_4.5_today,
D=dat_Nov$NN.sigma_today_2100_4.5,
E=dat_Nov$NN.sigma_2100_8.5_today,
F=dat_Nov$NN.sigma_today_2100_8.5)
write.csv(for_liqing, paste0(results_dir,"SigmaD_Hem_plot1.csv"), row.names=FALSE)
for_liqing_Md <- data.frame(No=dat_Nov$No, lat=dat_Nov$lat,long=dat_Nov$long,A=dat_Nov$NN.Mdist_today_1800, B=dat_Nov$NN.Mdist_1800_today, C=dat_Nov$NN.Mdist_2100_4.5_today, D=dat_Nov$NN.Mdist_today_2100_4.5, E=dat_Nov$NN.Mdist_2100_8.5_today, F=dat_Nov$NN.Mdist_today_2100_8.5)
write.csv(for_liqing_Md, paste0(results_dir,"SigmaD_Hem_plot2.csv"), row.names=FALSE)#plot(final_dat12$lat_today_2100_8.5, final_dat12$lat) #If A is today and B is future,
# where station No's climate in the future will come from today (or the closest similar climate today).
# Latitude of the station today (x-axis) where the latitude of the station in the future (y-axis) will come from
#plot(final_dat12$lat, final_dat12$lat_2100_8.5_today) # If A is future and B is today, this
# represents where station "No" will be found in the future (or the closest similar climate).
# Latitude of the station today (x-axis) and where it will be in the future (y-axis)Calculate percentages
### Calculate percentages ####
N <- nrow(dat_Nov)
calc_perc <- function(x, descrip){
print(paste("Percent of cells with moderate dissimilarity for", descrip,":",
round(sum(x>2 & x<4)/N,3)*100, sep=" "))
print(paste("Percent of cells with extreme dissimilarity for", descrip,":",
round(sum(x>4)/N, 3)*100, sep=" "))
}
calc_perc(dat_Nov$NN.sigma_today_1800, "1800 disappearing in 2000")## [1] "Percent of cells with moderate dissimilarity for 1800 disappearing in 2000 : 12.4"
## [1] "Percent of cells with extreme dissimilarity for 1800 disappearing in 2000 : 0"calc_perc(dat_Nov$NN.sigma_1800_today, "Novel climates in 2000 since 1800")## [1] "Percent of cells with moderate dissimilarity for Novel climates in 2000 since 1800 : 3.7"
## [1] "Percent of cells with extreme dissimilarity for Novel climates in 2000 since 1800 : 0"calc_perc(dat_Nov$NN.sigma_2100_4.5_today, "2000 disappearing in 2100 4.5")## [1] "Percent of cells with moderate dissimilarity for 2000 disappearing in 2100 4.5 : 34.2"
## [1] "Percent of cells with extreme dissimilarity for 2000 disappearing in 2100 4.5 : 35.6"calc_perc(dat_Nov$NN.sigma_today_2100_4.5, "Novel climates in 2100 4.5 since 2000")## [1] "Percent of cells with moderate dissimilarity for Novel climates in 2100 4.5 since 2000 : 42.1"
## [1] "Percent of cells with extreme dissimilarity for Novel climates in 2100 4.5 since 2000 : 10.3"calc_perc(dat_Nov$NN.sigma_2100_8.5_today, "2000 disappearing in 2100 8.5")## [1] "Percent of cells with moderate dissimilarity for 2000 disappearing in 2100 8.5 : 3.5"
## [1] "Percent of cells with extreme dissimilarity for 2000 disappearing in 2100 8.5 : 95"calc_perc(dat_Nov$NN.sigma_today_2100_8.5, "Novel climates in 2100 8.5 since 2000")## [1] "Percent of cells with moderate dissimilarity for Novel climates in 2100 8.5 since 2000 : 6.7"
## [1] "Percent of cells with extreme dissimilarity for Novel climates in 2100 8.5 since 2000 : 81.9"# Visualize relationship between sigma_D and M_D
plot(dat_Nov$NN.sigma_today_1800, dat_Nov$NN.Mdist_today_1800)
plot(dat_Nov$NN.sigma_today_2100_8.5, dat_Nov$NN.Mdist_today_2100_8.5, col=dat_Nov$numPCs_today_2100_8.5)
legend(0, 15, 1:4, 1:4)
abline(h=12)
hist(dat_Nov$numPCs_today_2100_8.5)